% M-File grxxm113.m
% Schiess-Verfahren fuer das lineare RWP
% u'' - exp(-t) u' -u =1, u(0) = 1, u(1) =1
% Loesung:
% u(t) = c1*(1-exp(t))+c2*exp(t-exp(-t))-exp(t)

clear;

% Grafikfenster (nur MATLAB)
screen=get(0,'ScreenSize');
screenplot=[50 50 screen(3)*0.8 screen(4)*0.8];
figure('Position',screenplot, ...
          'color',[1 1 0.77]);

clf;

% Schliessen Grafikfenster (nur MATLAB) 
schliessx=(screenplot(3)-screenplot(1))/2;
schliess=[schliessx-60 10 120 30];
uicontrol(gcf,...
         'Style','Push',... 
         'Position',schliess,...
         'String','Schliessen',...
         'Backgroundcolor',[0.5 0.93 0.9],...
         'callback','close');

% folgende Anweisungen arbeiten in MATLAB und OCTAVE:
% Exakte Loesung des RWP
tt=[0:1/100:1];
c1=(1+exp(1)-2*exp(1)^2/exp(1)^(1/exp(1)))/(1-exp(1));
c2=2*exp(1);
uu=c1*(1-exp(tt))+c2*exp(tt-exp(-tt))-exp(tt);
subplot(2,2,1);
plot(tt,uu,'b');
title('Loesung des RWP','FontSize',12);

% Lineares Schiess-Verfahren (R-K-verfahren) mit h=0.1
h=0.1;
y1(1,1:3)=[0,1,0];
for k=1:10
    y1(k+1,1)=k*h;
    k1=[y1(k,3);exp(-y1(k,1))*y1(k,3)+y1(k,2)+1];
    k2=[y1(k,3)+k1(2)*h/2;exp(-y1(k,1)-h/2)*(y1(k,3)+k1(2)*h/2)+y1(k,2)+k1(1)*h/2+1];
    k3=[y1(k,3)+k2(2)*h/2;exp(-y1(k,1)-h/2)*(y1(k,3)+k2(2)*h/2)+y1(k,2)+k2(1)*h/2+1];
    k4=[y1(k,3)+k3(2)*h;exp(-y1(k,1)-h)*(y1(k,3)+k3(2)*h)+y1(k,2)+k3(1)*h+1];
    w=[k1+2*k2+2*k3+k4]*h/6;
    y1(k+1,2)=y1(k,2)+w(1);
    y1(k+1,3)=y1(k,3)+w(2);
end;  
w1(1,1:3)=[0,0,1];
for k=1:10
    w1(k+1,1)=k*h;
    k1=[w1(k,3);exp(-w1(k,1))*w1(k,3)+w1(k,2)];
    k2=[w1(k,3)+k1(2)*h/2;exp(-w1(k,1)-h/2)*(w1(k,3)+k1(2)*h/2)+w1(k,2)+k1(1)*h/2];
    k3=[w1(k,3)+k2(2)*h/2;exp(-w1(k,1)-h/2)*(w1(k,3)+k2(2)*h/2)+w1(k,2)+k2(1)*h/2];
    k4=[w1(k,3)+k3(2)*h;exp(-w1(k,1)-h)*(w1(k,3)+k3(2)*h)+w1(k,2)+k3(1)*h];
    w=[k1+2*k2+2*k3+k4]*h/6;
    w1(k+1,2)=w1(k,2)+w(1);
    w1(k+1,3)=w1(k,3)+w(2);
end; 

% Lineare Kombination der zwei Loesungen
C1=(1-y1(11,2))/w1(11,2);
yy1=y1(:,2)+C1*w1(:,2);

subplot(2,2,2);
plot(tt,uu,'b',y1(:,1),yy1,'r');
title('Lineares Schiess-Verfahren, h = 0.1','FontSize',12);
legend('Loesung des RWP','Schiess-Verfahren, h=0.1');

% Lineares Schiess-Verfahren (R-K-verfahren) mit h=0.01
h=0.01;
y2(1,1:3)=[0,1,0];
for k=1:100
    y2(k+1,1)=k*h;
    k1=[y2(k,3);exp(-y2(k,1))*y2(k,3)+y2(k,2)+1];
    k2=[y2(k,3)+k1(2)*h/2;exp(-y2(k,1)-h/2)*(y2(k,3)+k1(2)*h/2)+y2(k,2)+k1(1)*h/2+1];
    k3=[y2(k,3)+k2(2)*h/2;exp(-y2(k,1)-h/2)*(y2(k,3)+k2(2)*h/2)+y2(k,2)+k2(1)*h/2+1];
    k4=[y2(k,3)+k3(2)*h;exp(-y2(k,1)-h)*(y2(k,3)+k3(2)*h)+y2(k,2)+k3(1)*h+1];
    w=[k1+2*k2+2*k3+k4]*h/6;
    y2(k+1,2)=y2(k,2)+w(1);
    y2(k+1,3)=y2(k,3)+w(2);
end;  
w2(1,1:3)=[0,0,1];
for k=1:100
    w2(k+1,1)=k*h;
    k1=[w2(k,3);exp(-w2(k,1))*w2(k,3)+w2(k,2)];
    k2=[w2(k,3)+k1(2)*h/2;exp(-w2(k,1)-h/2)*(w2(k,3)+k1(2)*h/2)+w2(k,2)+k1(1)*h/2];
    k3=[w2(k,3)+k2(2)*h/2;exp(-w2(k,1)-h/2)*(w2(k,3)+k2(2)*h/2)+w2(k,2)+k2(1)*h/2];
    k4=[w2(k,3)+k3(2)*h;exp(-w2(k,1)-h)*(w2(k,3)+k3(2)*h)+w2(k,2)+k3(1)*h];
    w=[k1+2*k2+2*k3+k4]*h/6;
    w2(k+1,2)=w2(k,2)+w(1);
    w2(k+1,3)=w2(k,3)+w(2);
end; 

% Lineare Kombination der zwei Loesungen
C2=(1-y2(101,2))/w2(101,2);
yy2=y2(:,2)+C2*w2(:,2);

subplot(2,2,3);
plot(tt,uu,'b',y2(:,1),yy2,'r');
title('Lineares Schiess-Verfahren, h = 0.01','FontSize',12);
legend('Loesung des RWP','Schiess-Verfahren, h=0.01');


% Lineares Schiess-Verfahren (R-K-verfahren) mit h=0.001
h=0.001;
y3(1,1:3)=[0,1,0];
for k=1:1000
    y3(k+1,1)=k*h;
    k1=[y3(k,3);exp(-y3(k,1))*y3(k,3)+y3(k,2)+1];
    k2=[y3(k,3)+k1(2)*h/2;exp(-y3(k,1)-h/2)*(y3(k,3)+k1(2)*h/2)+y3(k,2)+k1(1)*h/2+1];
    k3=[y3(k,3)+k2(2)*h/2;exp(-y3(k,1)-h/2)*(y3(k,3)+k2(2)*h/2)+y3(k,2)+k3(1)*h/2+1];
    k4=[y3(k,3)+k3(2)*h;exp(-y3(k,1)-h)*(y3(k,3)+k3(2)*h)+y3(k,2)+k3(1)*h+1];
    w=[k1+2*k2+2*k3+k4]*h/6;
    y3(k+1,2)=y3(k,2)+w(1);
    y3(k+1,3)=y3(k,3)+w(2);
end;  
w3(1,1:3)=[0,0,1];
for k=1:1000
    w3(k+1,1)=k*h;
    k1=[w3(k,3);exp(-w3(k,1))*w3(k,3)+w3(k,2)];
    k2=[w3(k,3)+k1(2)*h/2;exp(-w3(k,1)-h/2)*(w3(k,3)+k1(2)*h/2)+w3(k,2)+k1(1)*h/2];
    k3=[w3(k,3)+k2(2)*h/2;exp(-w3(k,1)-h/2)*(w3(k,3)+k2(2)*h/2)+w3(k,2)+k2(1)*h/2];
    k4=[w3(k,3)+k3(2)*h;exp(-w3(k,1)-h)*(w3(k,3)+k3(2)*h)+w3(k,2)+k3(1)*h];
    w=[k1+2*k2+2*k3+k4]*h/6;
    w3(k+1,2)=w3(k,2)+w(1);
    w3(k+1,3)=w3(k,3)+w(2);
end; 

% Lineare Kombination der zwei Loesungen
C3=(1-y3(1001,2))/w3(1001,2);
yy3=y3(:,2)+C3*w3(:,2);

subplot(2,2,4);
plot(tt,uu,'b',y3(:,1),yy3,'r');
title('Lineares Schiess-Verfahren, h = 0.001','FontSize',12);
legend('Loesung des RWP','Schiess-Verfahren, h=0.001');


% Tabelle 
disp('----------------------------------------------------------------------');
disp('               Lineares Schiess-Verfahren');
disp('----------------------------------------------------------------------');
disp(' t_n        exakt         h = 0.1         h = 0.01        h = 0.001     ');
disp('----------------------------------------------------------------------');

for s=0:10
    fprintf('%1.2f %15.10f %15.10f %15.10f %15.10f',...
        tt(1+s*10),uu(1+s*10),yy1(1+s),yy2(1+10*s),yy3(1+100*s));
    fprintf('\n');
end;
