function [y]=trapets(A,x,y0)
% Beräkna en approximation av lösningen till begynnelse
% värdes-problemet y'=Ay,  y(x(1))=y0  med hjälp av 
% trapetsmetoden och steglängden h=x(n+1)-x(n).
% A är en kvadratisk matris
%
N=length(x); h=x(2)-x(1);
y(:,1)=y0;
I=eye(size(A));     % enhetsmatrisen
[L,R]=lu(I-h/2*A);
for n=1:N-1         % Lös (I-h/2*A)y(n+1)=(I+h/2*A)y(n)
  y(:,n+1)=R\(L\(I+h/2*A))*y(:,n);
end

