clear all;
clc;

% Rozmiar siatki
K = input('K ');
L = input('L ');
%L = 100;
%K = 100;
h = 1/(K-1);
g = 1/(L-1);

% Punkty siatki
x = 0:g:1;
t = 0:h:1;

% Warunek początkowy
w_pocz = 0;
for j=1:L,
    y(1,j) = w_pocz;
end
for i=1:K,
    y(i,L) = 0;
end


% Rozwiązanie dokładne i ksi
for j=1:L
    for i=1:K
        yd(i,j) = sin(pi*x(j))*exp(-pi*pi*t(i))*(1-cos(4*pi*t(i)));
        ksi(i,j) = 4*pi*sin(pi*x(j))*sin(4*pi*t(i))*exp(-pi*pi*t(i));
    end
end

% Parametry pomocnicze
del = 1/2 - (g*g)/(12*h);
A = -del/(g*g);
C = A;
B = (2*del)/g^2 + 1/h;
alfa(2) = 0;
beta(2) = 0;

% Poszukiwanie rozwiązania
for i=1:K-1,
    for j=2:L-1,
            F(i,j) = -ksi(i,j)-((1-del)/g^2)*(y(i,j+1)+y(i,j-1))+((2-2*del)/g^2-1/h)*y(i,j);
            alfa(j+1) = C/(-B-alfa(j)*A);
            beta(j+1) = ((A*beta(j))+F(i,j))/(-B-alfa(j)*A);
    end
        for j=L-1:-1:1,
           y(i+1,j) = alfa(j+1)*y(i+1,j+1)+beta(j+1);
        end
end

% Błąd średniokwadratowy
suma = 0;
for i=1:K
    for j=1:L
        blad_q(i,j) = (yd(i,j)-y(i,j))^2;
        suma = suma+blad_q(i,j);
    end
end

blad = sqrt(suma)/(K*L);
mesh(x,t,y)