another function for solving differential equation other than dsolve in MATLAB

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])

[Rik]

(removed from spam filter)

[Bjorn Gustavsson]

You might get more and better response if you also write the ODEs in equation-form, the editing allows Latex-style equations, so you should be able to.

[Leila Dehdari]

Hi Bjorn, thanks for the response. Can you give me an example of your suggestion? if you mean diff(T,y) after dsolve function, this is the way that MATLAB help document shows.

[Bjorn Gustavsson]

No, I mean write the equations like you would see them in a book or article:

That makes it far easier for others to read than your long lines of code.

[Leila Dehdari]

Thanks Bjorn, I added it to the question

Answer 1

[Kiran Felix Robert]

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

[Leila Dehdari]

Thanks Kiran