function [x,y,t,u,A,f] = HeatEq
% Heat diffusion in a slab

% Warning: this code assumes case dependet os.

% Setup
L = 1.0;
W = 1.0;
T = 75;

Kx = 0.001;
Ky = Kx;

m = 10;
n = 10;
s = 75;

dx = L/m;
dy = W/n;
dt = T/s;

u = zeros((m+1)*(n+1),s+1);

A = zeros((m+1)*(n+1),(m+1)*(n+1));

f = zeros((m+1)*(n+1),1);

% Create mesh
t = (0:s)*dt;
x = (0:m)*dx;
y = (0:n)*dy;


alphax = Kx*dt/dx^2;
alphay = Ky*dt/dy^2;

alpha = 0;
beta = 1;
gamma = 0;


for r=1:n
    % Time stepping

    for j=2:n
        for i=2:m
            k = (j-1)*(m+1) + i;
            l = k;

            A(k,l-(m+1)) = alphay;
            A(k,l-1) = alphax;
            A(k,l) = 1 -(2*alphax + 2*alphay);
            A(k,l+1) = alphax;
            A(k,l+(m+1)) = alphay;
        end
    end

    j=1;
    i=1;
    k = (j-1)*(m+1) + i;
    l = k;

    A(k,l) = -3*alphay*dy/2*alpha/(-sqrt(2)) - 3*alphax*dx/2*alpha/(-sqrt(2)) + beta;
    A(k,l+(m+1)) = 4*alphay*dy/2*alpha/(-sqrt(2));
    A(k,l+2*(m+1)) = -alphay*dy/2*alpha/(-sqrt(2));
    A(k,l+1) = 4*alphax*dx/2*alpha/(-sqrt(2));
    A(k,l+2) = -alphax*dx/2*alpha/(-sqrt(2));
    
    f(k) = 1;

    for i=2:m
        k = (j-1)*(m+1) + i;
        l = k;

        A(k,l) = -3*alphay*dy/2*alpha + beta;
        A(k,l+(m+1)) = 4*alphay*dy/2*alpha;
        A(k,l+2*(m+1)) = -alphay*dy/2*alpha;

            f(k) = 1;

    end

    i=m+1;

    k = (j-1)*(m+1) + i;
    l = k;

    A(k,l) = -3*alphay*dy/2*alpha/(-sqrt(2)) + 3*alphax*dx/2*alpha/(sqrt(2)) + beta;
    A(k,l+(m+1)) = 4*alphay*dy/2*alpha/(-sqrt(2));
    A(k,l+2*(m+1)) = -alphay*dy/2*alpha/(-sqrt(2));
    A(k,l-1) = -4*alphax*dx/2*alpha/(sqrt(2));
    A(k,l-2) = alphax*dx/2*alpha/(sqrt(2));
    
        f(k) = 1;


    for j=2:n
        i = 1;
        k = (j-1)*(m+1) + i;
        l = k;
        A(k,l) = - 3*alphax*dx/2*alpha + beta;
        A(k,l+1) = 4*alphax*dx/2*alpha;
        A(k,l+2) = -alphax*dx/2*alpha;

            f(k) = 1;

        i=m+1;
        k = (j-1)*(m+1) + i;
        l = k;
        A(k,l) = 3*alphax*dx/2*alpha + beta;
        A(k,l-1) = -4*alphax*dx/2*alpha;
        A(k,l-2) = alphax*dx/2*alpha;
        
            f(k) = 1;

    end


    j=n+1;
    i=1;

    k = (j-1)*(m+1) + i;
    l = k;

    A(k,l) = 3*alphay*dy/2*alpha/(sqrt(2)) - 3*alphax*dx/2*alpha/(-sqrt(2)) + beta;
    A(k,l-(m+1)) = -4*alphay*dy/2*alpha/(sqrt(2));
    A(k,l-2*(m+1)) = alphay*dy/2*alpha/(sqrt(2));
    A(k,l+1) = 4*alphax*dx/2*alpha/(-sqrt(2));
    A(k,l+2) = -alphax*dx/2*alpha/(-sqrt(2));

        f(k) = 1;

    for i=2:m
        k = (j-1)*(m+1) + i;
        l = k;

        A(k,l) = 3*alphay*dy/2*alpha + beta;
        A(k,l-(m+1)) = -4*alphay*dy/2*alpha;
        A(k,l-2*(m+1)) = alphay*dy/2*alpha;
        
            f(k) = 1;

    end

    i=m+1;

    k = (j-1)*(m+1) + i;
    l = k;

    A(k,l) = 3*alphay*dy/2*alpha/(sqrt(2)) + 3*alphax*dx/2*alpha/(sqrt(2)) + beta;
    A(k,l-(m+1)) = -4*alphay*dy/2*alpha/(sqrt(2));
    A(k,l-2*(m+1)) = alphay*dy/2*alpha/(sqrt(2));
    A(k,l-1) = -4*alphax*dx/2*alpha/(sqrt(2));
    A(k,l-2) = alphax*dx/2*alpha/(sqrt(2));
    
        f(k) = 1;


end