function [t,U]=theta_method1(lambda,N,theta)

% THETA_METHOD1 solve ODE with theta method
%   [t,U]=theta_method1(lambda,N,theta) returns a sequence of time points t
%   (in the interval [0,1]) and an approximate solution U at those time 
%   points to the problem u'(t)=lambda*u with u(0)=1 (with exact solution
%   u(t)=exp(lambda*t)).
%   The approximate solution is computed using the theta method with N+1
%   timesteps.

t=linspace(0,1,N+1);
dt=t(2)-t(1);

U=zeros(1,N+1);
U(1)=1;

for n=1:N
    U(n+1)=U(n)*(1+(1-theta)*lambda*dt)/(1-theta*lambda*dt);
end