MATLAB: Another function for solving differential equation other than dsolve in MATLAB

dsolve

I am using dsolve function in MATLAB to solve a differential equation. It is so slow and can not end up with a result. Is there a way to speed it up, or is there another way for solving a differential equation? I don't want to solve it numerically.
The equations are as follow:
all the derivatives in above equations are replaced with equivalent expressions, so all are only function of y and T only.
I need to solve this equation:
Cg, Cs and qi are constants. The initial conditions for positive root is (y=0.01, T=0) and for negative root is (y=0.99, T=0).
Here is my code:
R=1.9872;                          %cal/mol/K
eps=0.763;
rhos=0.484;                        %gr/cm3
T0=300;                            %theta K
y0=0.01;
ys=0.99;
p=1;                              %atm
Upg=0.68e-3;                      %mol/cm2/s
L=100;                            %cm   column length
Cg=1.04*239.006/1000*28.0134;      % kj/kg/k  to  cal/mol/k   N2
Cs=0.22;                          %cal/gr/k
%  m(mol/gr)     b(1/atm)        q(cal/mol)
cA=[3.65e-3      2.8e-4            4900];  %  co2 BPL AC
cB=[3.65e-3      13.5e-4           2500];  %  N2 BPL AC 
i=1;
IPA=cA(i,:); IPB=cB(i,:);
IPA1=IPA(1);IPA2=IPA(2);IPA3=IPA(3);
IPB1=IPB(1);IPB2=IPB(2);IPB3=IPB(3);
syms T(y)
DN1T= @(y,T)    (IPA1*IPA2*p*y*exp(IPA3/(R*(T + T0)))*((IPA2*IPA3*p*y*exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2) - (IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))...
    *(y - 1))/(R*(T + T0)^2)))/(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2 -...
    (IPA1*IPA2*IPA3*p*y*exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2*(IPA2*p*y*exp(IPA3/(R*(T + T0))) - ...
    IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1));
DN2T =@(y,T)     (IPB1*IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))*(y - 1))/(R*(T + T0)^2*(IPA2*p*y*exp(IPA3/(R*(T + T0))) - ...
    IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)) - (IPB1*IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1)*((IPA2*IPA3*p*y*...
    exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2) - (IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))*(y - 1))/(R*(T + T0)^2)))...
    /(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2;
 
 
DN1Y =@(y,T)      (IPA1*IPA2*p*exp(IPA3/(R*(T + T0))))/(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)...
    - (IPA1*IPA2*p*y*exp(IPA3/(R*(T + T0)))*(IPA2*p*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))))/...
    (IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2;
 
DN2Y =@(y,T)      (IPB1*IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1)*(IPA2*p*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))))/...
    (IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2 - (IPB1*IPB2*p*exp(IPB3/(R*(T + T0))))...
    /(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1);
  
A=@(y,T)   Cg*(DN1T(y,T)-y*(DN1T(y,T)+DN2T(y,T)));
B=@(y,T)   (Cg*(T)*(DN1T(y,T)+DN2T(y,T))-Cs)+Cg*(DN1Y(y,T)-y*(DN1Y(y,T)+DN2Y(y,T)))+cA(i,3)*DN1T(y,T)+cB(i,3)*DN2T(y,T);
C=@(y,T)    Cg*(T)*(DN1Y(y,T)+DN2Y(y,T))+cA(i,3)*DN1Y(y,T)+cB(i,3)*DN2Y(y,T);
 
 S4=@(y,T)   (-B(y,T) +(B(y,T) ^2-4*A(y,T) *C(y,T) )^0.5)/2/A(y,T) ;
 S2=@(y,T)   (-B(y,T) -(B(y,T) ^2-4*A(y,T) *C(y,T) )^0.5)/2/A(y,T) ;
S4sol=dsolve([diff(T,y)==S4(y,T) , T(y0)==0])
S2sol=dsolve([diff(T,y)==S2(y,T) , T(ys)==0])

Best Answer

Hi Leila,
One way to speed-up the execution is to add a limit on the maximum degree of radicals. The following shows you an example, 
S = dsolve(..,'MaxDegree',2);
This will force dsolve to assume implicit formulas for polynomials of degree greater than the specified value.
Refer to the dsolve documentation for more information.
Kiran Felix Robert
Related Question