MATLAB: Problem solving a system of differential equations

differential equationsdsolveMATLABode

I am trying to solve a system of differential equations in Matlab.
   dn/du=(-2*u*n-K*(n*u-(1+g)))/(1+u^2+K*u*(u-(1+g)/n))
   dxi/du=(1-u^2)/(1+u^2+K*u*(u-(1+g)/n))
   df/du=(2*u+K*u^2*(u-(1+g)/n))/(1+u^2+K*u*(u-(1+g)/n))
K and gamma are constants. and i have the initial conditions ( n(0)=1, xi(0)=0, f(0)=0 )
First i tried to use 'dsolve' as shown below.
 g=0.1;
 K=3;
 syms n(u) u 
 n(u)=dsolve(diff(n,u)== (-2*u*n-K*(n*u-(1+g)))/(1+u^2+K*u*(u-(1+g)/n)),n(0)==1)
 syms x(u) u n
 x(u)= dsolve(diff(x,u)== (1-u^2)/(1+u^2+K*u*(u-(1+g)/n)),x(0)==0)
 syms f(u) u n
 f(u)=dsolve(diff(f,u)== (2*u+K*u^2*(u-(1+g)/n))/(1+u^2+K*u*(u-(1+g)/n)),f(0)==0)
from this i got that there is not explicit solution for the first equation and returned an [empty sym] as shown below.
 Warning: Explicit solution could not be found.
 In dsolve (line 201)
 In Untitled (line 37)
 n(u) = [ empty sym ]
 x(u) = (1089*atan(33/(1600*n^2 - 1089)^(1/2) - (80*n*u)/(1600*n^2 - 1089)^(1/2)))/(160*n*(1600*n^2 - 1089)^(1/2)) - (35937*log(40*n*u^2 - 33*u + 10*n))/(2*(- 256000*n^3 + 174240*n)) - u/4 + (33*log(10*n)*(1600*n^2 - 1089)^(1/2) - 2178*atan(33/(1600*n^2 - 1089)^(1/2)) + 8000*n^2*atan(33/(1600*n^2 - 1089)^(1/2)))/(320*n*(1600*n^2 - 1089)^(1/2)) - (25*n*atan(33/(1600*n^2 - 1089)^(1/2) - (80*n*u)/(1600*n^2 - 1089)^(1/2)))/(1600*n^2 - 1089)^(1/2) + (26400*n^2*log(40*n*u^2 - 33*u + 10*n))/(- 256000*n^3 + 174240*n)
 f(u) = (3*u^2)/8 - (33*u)/(160*n) - (log(40*n*u^2 - 33*u + 10*n)*(3200000*n^4 - 3920400*n^2 + 1185921))/(2*(-    10240000*n^4 + 6969600*n^2)) + (log(10*n)*(3200000*n^4 - 3920400*n^2 + 1185921))/(2*(- 10240000*n^4 + 6969600*n^2)) +  (33*atan((1089*(2800*n^2 - 1089))/((1600*n^2 - 1089)^(1/2)*(92400*n^2 - 35937)))*(2800*n^2 - 1089))/(6400*n^2* (1600*n^2 - 1089)^(1/2)) - (33*atan((1089*(2800*n^2 - 1089))/((1600*n^2 - 1089)^(1/2)*(92400*n^2 - 35937)) - (2640*nu(2800*n^2 - 1089))/((1600*n^2 - 1089)^(1/2)*(92400*n^2 - 35937)))*(2800*n^2 - 1089))/(6400*n^2*(1600*n^2 - 1089)^(1/2))
Then i tried to use 'ode45'.
 u=0:0.1:1;
 K=3;
 g=0.1;
 dndu=@(u,n) (-2*u*n-K*(n*u-(1+g)))/(1+u^2+K*u*(u-(1+g)/n));
 [u,n]=ode45(dndu, u, 1) % initial n=1
 dxidu=@(u,xi) (1-u^2)/(1+u^2+K*u*(u-(1+g)/n));
 [u,xi]=ode45(dxidu, u, 0); %initial xi=0

 dfdu=@(u,f) (2*u+K*u^2*(u-(1+g)/n))/(1+u^2+K*u*(u-(1+g)/n));
 [u,f]=ode45(dfdu, u, 0) %initial f=0
K and g are constants and u gets values for 0 to 1. When i run this the first equation is solved and got an answer for u and n but i get an error for the second equation that matrix dimension must agree as shown below.
 u = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
 n = 1.0000 1.3409 1.6243 1.7746 1.7945 1.7263 1.6119 1.4800 1.3470 1.2208 1.1048
 Error using /
 Matrix dimensions must agree.
 Error in @(u,xi)(1-u^2)/(1+u^2+Ku(u-(1+g)/n))
 Error in odearguments (line 87)
 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
 Error in ode45 (line 113)
 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
 Error in ode45method (line 24)
 [u,xi]=ode45(dxidu, u, 0); %initial xi=0
Can anyone help me how to approach this type of problem? Thanks in advance!

Best Answer

uspan=0:0.1:1;
K=3;
g=0.1;
dydu=@(u,y)[(-2*u*y(1)-K*(y(1)*u-(1+g)))/(1+u^2+K*u*(u-(1+g)/y(1)))
            (1-u^2)/(1+u^2+K*u*(u-(1+g)/y(1)))
            (2*u+K*u^2*(u-(1+g)/y(1)))/(1+u^2+K*u*(u-(1+g)/y(1)))];
[U,Y]=ode45(dydu,uspan,[1 0 0]);
plot(U,Y(:,1),U,Y(:,2),U,Y(:,3));
Best wishes
Torsten.
Related Question