Suppose that your list of equations is copied from Matlab as a single string, which I'll call MS (Matlab String), like this

MS:= "C_p_e = C_state/C_c;
C_p_f = I_state/I_i;
R_p_e = R_r*C_p_f;
I_p_e = (Se_p_e-R_p_e)-C_p_e;";

The presence of newline characters in the string is insigificant---Maple ignores them. The semicolons, however, are necessary. Then do

eq:= parse~(StringTools:-Split(MS, ";"));

Now eq is a list of all the equations, so eq[1] is C_p_e = C_state/C_c, etc. You just need to use subs once, on all of eq. It's more efficient than calling allsubs repeatedly.

eq_a:= subs(
     C_p_e=fs(t), C_p_f=v(t), R_p_e=fd(t), I_p_e=fm(t), C_c=1/k, I_i=m, 
     R_r=c, Se_p_e=fe(t), I_state=int(fm(t),t), C_state=int(v(t),t),
     eq
);

The solutions are not necessarily real in your example. If the solutions were necessarily real, you could use syntax almost identical to the Mathematica:

a:= INTERVAL(2..4):
b:= INTERVAL(10..11):
c:= INTERVAL(1..3):
evalr~([solve(a*x^2+b*x+c)]);

If you want to do complex range arithmetic with your original example, you could do this:

a:= INTERVAL(2,0,4,0):
b:= INTERVAL(1,0,3,0):
c:= INTERVAL(1,0,3,0):
evalrC~([solve(a*x^2+b*x+c)]);

However, the results are not perfect, because they are in floating point. Note how a complex interval is specified with four numbers. See ?evalrC and ?evalr.

simplify(x1=a-y1-d*y2, {a-y2-d*y1= x2, 1-d^2= b, a-a*d= c});

Note that the three substitution equations must be written with the expressions on the left and the single variables on the right.

Don't use map when applying expand or simplify to a single expression. Use map to distribute these operations over a Matrix, but not on an single entry extracted from a Matrix.

There are also "tutors" for this in the Student package, but I think that it's instructive to see the components of a plot command.

Download RatFuncPlot.mw

In the future, you should show what you tried. Here's how to plot it.

f:= x-> 2*cos(x)-x^2:  a:= 2:  r:= 3:
plot([f(x), D(f)(a)*(x-a)+f(a)], x= a-r..a+r, color= [red,blue], linestyle= [solid,dot]);

I hope that you can generalize it from there.

Here's a solution for an arbitrary number of dancers, using evenly spaced points on the unit circle as the starting points. The ODEs are essentially the same as in my solution above. I decided to form the dependent variable names with indexing (x[1], x[2], y[1], etc.) rather than with cat and ||. I also changed the indexing from 0..n-1 to 1..n. The brevity of the code shows Maple's great power at forming and working with a repetitive set of ODEs.

restart:
SpiralDance:= proc(n::posint)
local
     k, K:= k= 1..n, x, y, V_k:= [x[k],y[k]], v,
     eqns:= seq(((D+(X->X))(v[k]) = v[1+irem(k,n)]) $ K, v= [x,y]),
     ics:= op~([(V_k =~ [cos,sin](-2*k*Pi/n)) $ K])[],
     t, Sol:= dsolve({eqns(t),ics(0)}, numeric)
;
     plots:-odeplot(Sol, [V_k(t) $ K], thickness= 3, axes= none, gridlines= false)
end proc:

SpiralDance(5);

The breakpoint command in Maple is stopat. See ?stopat.

Your transform is (x,y)-> [x-y]. This says "take a point (x,y) on the graph and replace it with the single number x-y." That makes no sense: You need to replace points with points---pairs of numbers---not with single numbers.

I also encourage you to replace PLOT with plots:-display. 

There's a decimal point that you included somewhere in your input. Remove it.

Here's a simple solution, off the top of my head, with the dancers starting at [1,1], [1,-1], [-1,-1], [-1,1]. I leave animating this solution to you. It's really easy. The position of dancer k (0 to 3) at time t is [x||k(t), y||k(t)].

Download Spriral_Dance.mw

Try the old profiling commands: exprofile, excallgraph, profile. I don't have much experience with profile. I've used exprofile and excallgraph a lot, and I've never come across code that they couldn't handle.

The mathematical constant Pi is spelled with a capital P in Maple. If you spell it pi, it's just another variable, which, unfortunately, prettyprints exactly the same as the mathematical constant. However, the printed forms can be distinguished using the lprint command.

A (nonconstant) function of a random variable is itself a random variable. So your Y is already a random variable and there's no need to apply RandomVariable to it. Skip T1, and do what you want with Y, such as 

simplify(PDF(Y,t));

or, more directly,

simplify(PDF(1/X,t));

Download Pivot.mw