I have implemented following equations both with the Matlab own ODE solvers and with a very simple finite difference scheme. The latter works properly, while the ODE code does not produce suitable solutions. What I am doing wrong? These are a set of coupled implicit differential equations.
% Define the range of r for the solution
r_span = [r2, r_outlet];
% Initial state vector [cm, cu, rho, p]
y0 = [cm2, cu2, rho2, p2];
% Define the ODE system
syms cm(r) cu(r) rho(r) p(r)
eq2 = cm*r*diff(rho) + rho*r*diff(cm)+cm*rho == 0;
eq3 = diff(r*cu) == -r^2*cu*cf*rho*sqrt(cu^2+cm^2)/(rho2*cm2*b2*r2);
eq4 = diff(p)/rho == -cu*diff(cu)-cm*diff(cm)-(cd*rho*r* (cu^2+cm^2)^1.5)/(rho2*cm2*b2*r2);
eq5 = (gamma/(gamma-1))*diff(p/rho) == -cu*diff(cu)-cm*diff(cm);
% eqns = [eq1; eq2; eq3; eq4; eq5];
eqns = [ eq2; eq3; eq4; eq5];
[V, S] = odeToVectorField(eqns);
odefcn = matlabFunction(V, 'vars', {'r', 'Y'});
opt = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
% Solve the ODE system using ode45
[r_sol, y_sol] = ode15s(@(r, Y) odefcn(r, Y), r_span, y0,opt);
Thanks in advance!
As said before a comparison with a finite difference (FD) scheme, showed that these code cannot solve the system of equations. Also with FD it was possible to generate physically sound results. For example, conservation of mass is given for FD. In the ODE code this conservation is defined in eq1