function thetasystem(n,th)
% Solves the ivp y'=y +g , y_1=a, y_2=b 
% by the theta method
%
fprintf('\n')
%
for iit=1:10
%
nn=2^(iit-1)*n;
uu(iit)=1/nn;
%
h=1/nn;
%
A1=zeros(2,2);A2=zeros(2,2);y=zeros(2,1);er=zeros(nn+1,2);
ex=zeros(2,1);g1=zeros(2,1);g2=zeros(2,1);
%y(1,1)=1;y(2,1)=2;
y(1,1)=2;y(2,1)=0;
%
for i=0:nn
t=i*h;
t1=(i+1)*h;
ex=exact1(t);
er(i+1,:)=abs(ex(:,1)-y(:,1));
%
g1=gterm1(t);
g2=gterm1(t1);
A1=matrix1(t);
A2=matrix1(t1);
%
%y=(eye(2)+h*A)*y+h*g;
y=inv(eye(2)-h*(1-th)*A2)*(eye(2)+h*th*A1)*y+...
    h*inv(eye(2)-h*(1-th)*A2)*(th*g1+(1-th)*g2);
end
u1(iit)=max(er(:,1));
u2(iit)=max(er(:,2));
fprintf('h =%12.6f        maxerr1 = %15.6e    maxerr2 = %15.6e\n',h,u1(iit),u2(iit))
end
u1';
loglog(uu,u1,'b+',uu,u2,'ro');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exact=exact1(t);
exact=zeros(2,1);
%exact(1,1)=cos(t)+sin(t);
%exact(2,1)=2*cos(t);
exact(1,1)=exp(3*t)+exp(-t);
exact(2,1)=2*exp(3*t)-2*exp(-t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gterm=gterm1(t);
gterm=zeros(2,1);
%gterm(1,1)=4*cos(t)-2*sin(t);
%gterm(2,1)=5*cos(t)-5*sin(t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function matrix=matrix1(t);
matrix=zeros(2,2);
%matrix(1,1)=1;matrix(1,2)=-2;
%matrix(2,1)=3;matrix(2,2)=-4;
matrix(1,1)=1;matrix(1,2)=1;
matrix(2,1)=4;matrix(2,2)=1;