A data is not calculated, while using rk4 code

function [tp,yp] = eulersys(dydt, tspan,y0,h,varargin)
% rk4sys : fourth-order Runge-Kutta for a system of ODEs
%   [t,y] = rk4sys(dydt, tspan,y0,h,p1,p2,...) : integrates a system of
%   ODEs with fourth-order RK method
% input :
%   dydt = name of the M-file that evaluates the ODEs
%   tspan = [ti, tf]; initial and final times with output generated at
%                     interval of h, or
%   = [t0 t1 ... tf]; sepcific times wher solution output
%   y0 = initial values of dependent varaibles
%   h = step size
%   p1,p2,... = additiional parameters used by dydt
% output:
%   tp = vector of independent variable
%   yp = vector of solution for dependent variables

if nargin<4,error('at least 4 input arguments required'),end
ti = tspan(1);tf = tspan(2);
if ~(tf>ti),error('upper limit must be greater than lower'),end
t = (ti:h:tf)'; n = length(t);
% if necessary, add an additional value of t
% so that range goes from t = ti to tf
if t(n)<tf
t(n+1) = tf;
n = n+1;
end

tt = ti; y(1,:) = y0;
np = 1 ; tp(np) = tt ; yp(np,:) = y(1,:);
i = 1;
while(1)
    tend = t(np+1);
    hh = t(np+1) - t(np) ;
    if hh>h, hh = h ;end
    while(1)
        if tt+hh>tend, hh = tend-tt ; end
        k1 = dydt(tt,y(i,:), varargin{:})';
        y(i+1,:) = y(i,:) + k1.*hh;
       
        tt = tt+hh;
        i = i+1;
        if tt>=tend,break,end
    end
    
    np = np+1; t(np) = tt; yp(np,:) = y(i,:);
    if tt>=tf,break,end
end

i used the code above, which i rewrote from rk4sys, to solve an euler(rk1) question

but the result shows

>> [t y]=eulersys(@dydt,[0 0.4],[2 4],0.1)

t =

0

y =

2.0000 4.0000

2.1000 2.4000

2.1324 1.7952

2.1153 1.4516

2.0626 1.2287

the y datas are correct

but why is t not calculated/showed?

am i doing anything wrong with the function code?

Answer 1

[Adam Danz]

"why is t not calculated/showed"

It is shown. It equals 0. If you run the line below the extra space between the variable name and its value will be reduced.

format compact

The tp output will always be a scalar value (following your description of tspan).

Here's a reduced version of your code showing where the tp value comes from. tp will always equal the first element of tspan.

ti = tspan(1);
tt = ti;
np = 1 ; 
tp(np) = tt;

[Chm Lin]

thank you! i think i know how to do it:)