close all
clear all
clc
%liczba punktów na brzegu
n=18;
A=zeros(n^2,n^2);
b=ones(n^2,1);
h=1/(n+1); %hx=hy
%czas
c = 0.7;
dt = 0.01; %krok czasowy
tk=5; %czas końcowy -inaczej czas animacji
y=[];
%p - poprzedni, n - następny
up = zeros(n^2,1);
u = zeros(n^2,1);
un = zeros(n^2,1);
figure;

%macierz A

for t=0:dt:tk
    %dla punktów na brzegu
    for i=1:n^2
        if i<=n
            A(i,i)=1;
        elseif i>(n*n-n)
            A(i,i)=1;
        elseif mod(i,n) == 1
            A(i,i)=1;
        elseif mod(i,n) == 0
            A(i,i)=1;
        else
            %dla punktów wewnątrz obszaru
            A(i,i)=-4*c^2*dt^2-h^2;
            A(i,i-1)=c^2*dt^2;
            A(i,i+1)=c^2*dt^2;
            A(i,i-n)=c^2*dt^2;
            A(i,i+n)=c^2*dt^2;

        end
    end
%wektor b dla punktów na brzegu
    for i=1:n^2
        if i<=n
            b(i,1)=0;
        elseif i>(n*n-n)
            b(i,1)=0;
        elseif mod(i,n) == 0
            b(i,1)=0;
        elseif mod(i,n) == 1
            b(i,1)=0;
        else
            b(i,1)=-2*h^2*u(i)+up(i)*h^2;
        end
    end

    x1=[];
    y1=[];
    un=gmres(A,b,8,1e-9,100);
    y=reshape(u,n,n);
    surf(y);
    axis([0 n+1 0 n+1 -1 1]);
    shading interp
    colormap();
    pause(0.01);
    up=u;
    u=un;
    u(153,1)=sin(5*t);
end