MATLAB: ODE45 solver problem outcome

ode45odesolver

I have 5 differential equations and a couple of extra equations that I want to solve for 5 variables.
The problem is that if you look at the graph you see that the outcomes do not match. dh/dt should be -u1(t) but that is not the case in the graphs. I do not know where it goes wrong in my code.
If I put the equations for u1(t), u2(t) and ud(t) in the Eqns the same results are obtained which seems strange to me. How do I write the code such that all equations are true?
clear all;
ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13
syms u1(t) u2(t) ud(t) h(t) wr(t) T Y
u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == (u2(t)*wr(t))/wd;
Eqns = [diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(h(t),t) == - u1(t);
        diff(wr(t),t) == -us];
[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = linspace(0, 200, 251);
Y0 = [0; 0; 0; 0; 30]+0.001;
[T, Y] = ode45(DEfcn, tspan, Y0);

Best Answer

The ‘==’ are correct here. They are symbolic equations, not logical operations. Otherwise, I would agree.
Also, your code works correctly. The integrated value ‘h(t)’ will appear different from ‘-u1(t)’ because it is integrated. If you plot ‘u1(t)’ (integrated as ‘Y(:,2)’) against the negative of the derivative of ‘h(t)’: 
figure
plot(T,Y(:,2), '-b',  T,-gradient(Y(:,4),tspan(2)), '--r')
grid
(with the derivative calculated by the gradient function), they are almost exactly the same.
.
Related Question