the code in octave using theta method to solve an ode doesn't run ( can't find the error)

Using the theta method seen in the picture below

(https://i.sstatic.net/MJx0d.png)

for theta in [0.1].

I'm trying to solve some ODE seen in the code

function [tnodes,y] = theta_method(f, t_interval, N, y0, theta)

global glob_h glob_yold glob_tnew glob_f glob_theta;

y= y0;
yold= y0(:); %initial conditions
t0 = t_interval(1);
tN = t_interval(2);

tnodes = linspace(t0, tN, N+1); %the nodes of the partition

glob_h = (tN- t0)/N; %mesh-size
glob_f = f;
glob_theta=theta;
for tnew = tnodes(2:end)
  glob_tnew = tnew;
  glob_yold = yold;
  ynew = fsolve(@(x)F(x), yold);
  y =[y; ynew.'];
  yold = ynew;
end
endfunction

function rhs = F(x)

  global glob_h glob_yold glob_tnew glob_f glob_theta;

  rhs = x - glob_h*glob_theta*feval(glob_f, glob_tnew, x)-glob_h*(1.0 - glob_theta)*feval(glob_f, t0, glob_yold) - glob_yold;

Endfunction

And the plot :

PLOT
hold off
N=20;
t_interval = [0 1];
y0 = 1;
theta = 0.5;
[t,y]=theta_method(@model_f, t_interval, N, y0, theta)
uexact = model_exact(t)
plot(t,y,'linewidth', 1, 'g*'); % plot approx vs t
hold on
plot(t,uexact,'linewidth', 1, 'r'); % plot exact solution vs t
legend('BE approx','exact solution')
title('Solution of y''=\lambda y')
xlabel('time')

MODEL EXACT
function rhs = model_exact(t)
  lambda = -1;
  rhs = exp(lambda*t);
endfunction
MODEL F

function rhs = model_f(t,y)
  lambda = -1;
  rhs = lambda*y;
endfunction

the program says the the error is `

error: 't0' undefined near line 28, column 99
error: called from
    theta_method>F at line 28 column 7
    theta_method>@<anonymous> at line 18 column 21
    fsolve at line 246 column 8
    theta_method at line 18 column 8
    a at line 6 column 6

but i can't understand why what is wrong. I have define t0 and i think there isn't any syntax error.

Solution

It would be simpler to avoid the need for global variables.

told = tnodes(1)
for tnew = tnodes(2:end)
  fold = f(told,yold);
  yhalf = yold+(1-theta)*h*fold;
  F = @(x) yhalf+theta*h*f(tnew,x)-x;
  ynew = fsolve(F, yold+h*fold);
  y =[y; ynew.'];
  told = tnew;
  yold = ynew;
end

It would be advisable to scale the error tolerances of fsolve with h or h^2.

See SIR model using fsolve and Euler 3BDF for a similar construction for a different implicit method. Another method with setting the options for the error tolerances, Matlab code for $y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$