Errors when trying to plot a solution to a system of ODE

3 views (last 30 days)
Cris19 on 17 Dec 2019
Commented: Cris19 on 17 Dec 2019
I try to plot on [1000, 5000] the solution of the system of ODEs
with the initial conditions , where and .
I used the function
function dzdt=odefunw1(t,z)
f=1/(t+1);
g=1+exp(-t);
h=diff(f);
dzdt=zeros(2,1);
dzdt(1)=z(2)-f*z(1);
dzdt(2)=(g+h)*z(1)-f*z(2);
end
and the commands
tspan = [1000 5000];
z0 = [0.001 0.001];
[t,z] = ode45(@(t,z) odefunw1(t,z), tspan, z0);
plot(t,z(:,1),'r')
The following errors occured:
In an assignment A(:) = B, the number of elements in A and B must be the same.
Error in odefunw1 (line 7)
dzdt(2)=(g+h)*z(1)-f*z(2);
Error in @(t,z)odefunw1(t,z)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
How could I fix it ? Many thanks in advance.

Walter Roberson on 17 Dec 2019
[t,z] = ode45(@(t,z) odefunw1(t,z), tspan, z0);
Okay, ode45 will invoke odefunw1
function dzdt=odefun1(t,z)
which will receive the parameters through the same names as in the outside (no complications about different variable names.)
ode45 always passes the first parameter, t, as a numeric scalar.
f=1/(t+1);
Using the / operator with a scalar on the right hand side is acceptable and will produce the same value as if you had used the ./ operator, so this line is okay
g=1+exp(-t);
Another scalar, not a problem
h=diff(f);
diff() applied to a numeric array is the numeric difference function that calculates the numeric difference between adjacent entries. You are passing it a numeric scalar. Because there are no adjacent entries to a scalar, the output of diff() applied to a numeric scalar is empty.
dzdt(2)=(g+h)*z(1)-f*z(2);
Because h is empty, the right hand side of that equation is empty.
When you look back at the mathematical formula, we see that your h corresponds to the formula term where . You should be taking the mathematic derivative of that, getting so inside the function you should calculate h as that formula.
Cris19 on 17 Dec 2019
Thank you so very much for the thorough answer!