function ev = a1b(n,D,deltat,deltax)
    a = D*deltat/(deltax*deltax)
    for i=1:n
        for j=1:n
            if i == j then
                A(i,j) = 1-2*a
            else
                if abs(i-j) == 1 then
                    A(i,j) = a
                end
            end
        end
     end
     disp(A)
     ev = spec(A)
endfunction



function u = applyscheme(D,deltat,deltax,f,a1,aM,iter)
    n = length(f)
    u = 1:n
    
    for j=1:iter  
        u(1) = f(1) 
        for i=2:n-1
            u(i) = f(i)+D*deltat/(deltax*deltax)*(f(i+1)-2*f(i)+f(i-1))
        end
        u(n) = u(n-1)
        //u(n) = f(n)
        f = u
    end
endfunction

function a2a()
    clf()
    deltax = 0.05
    D = 1
    m = floor(0.5/deltax)+1
    disp(m)
    f(1) = 0
    f(m) = 1
    for i=2:m-1
        f(i)=sin(%pi*(i-1)*deltax)
    end
    
    deltat = 0.00125
    t = 0.07425
    iter = floor(t/deltat)
    u = applyscheme(D,deltat,deltax,f,0,0,iter)
    plot2d(f)
    plot2d(u)
    t = 0.111375
    u2 = applyscheme(D,deltat,deltax,f,0,0,floor(t/deltat))
    plot2d(u2)
    t = 0.136125
    u3 = applyscheme(D,deltat,deltax,f,0,0,floor(t/deltat))
    plot2d(u3)
    //disp(f)
    //disp(u)
    //disp(u2)
    //disp(u3)
endfunction


function a2b()
    clf()
    deltax = 0.05
    D = 1
    m = floor(0.5/deltax)+1
    disp(m)
    f(1) = 0
    f(m) = 1
    for i=2:m-1
        f(i)=sin(%pi*(i-1)*deltax)
    end
    
    deltat = 0.00132
    t = 0.07425
    iter = floor(t/deltat)
    u = applyscheme(D,deltat,deltax,f,0,0,iter)
    plot2d(f)
    plot2d(u)
    t = 0.111375
    u2 = applyscheme(D,deltat,deltax,f,0,0,floor(t/deltat))
    plot2d(u2)
    t = 0.136125
    u3 = applyscheme(D,deltat,deltax,f,0,0,floor(t/deltat))
    plot2d(u3)
    //disp(f)
    //disp(u)
    //disp(u2)
    //disp(u3)
endfunction

function u = cranknicolson(D,deltat,deltax,f,a1,aM,iter,omega)
    a = D*deltat/(deltax*deltax)
    disp(a)
    n = length(f)
    for j=1:iter 
    A = eye(n,n)
    for i=2:n-1
        for j=1:n
            if i == j then
                A(i,j) = 1-2*(1-omega)*a
            else
                if abs(i-j) == 1 then
                    A(i,j) = a*(1-omega)
                end
            end
        end
     end
     
     d = A*(f)
     A = eye(n,n)
    for i=2:n-1
        for j=1:n
            if i == j then
                A(i,j) = 1+2*omega*a
            else
                if abs(i-j) == 1 then
                    A(i,j) = -a*omega
                end
            end
        end
     end
    
     u = inv(A)*(d)
     
     f = u
     end
endfunction

function a3a(omega)
    clf()
    deltax = 0.05
    D = 1
    m = floor(1/deltax)+1
    disp(m)
    f(1) = 0
    f(m) = 0
    for i=2:m-1
        f(i)=sin(%pi*(i-1)*deltax)
    end
    
    deltat = 0.01
    
    disp(deltat)
    t = 0.07425
    iter = floor(t/deltat)
    u = cranknicolson(D,deltat,deltax,f,0,0,iter,omega)
    plot2d(f)
    plot2d(u)
    t = 0.111375
    u2 = cranknicolson(D,deltat,deltax,f,0,0,floor(t/deltat),omega)
    plot2d(u2)
    t = 0.136125
    u3 = cranknicolson(D,deltat,deltax,f,0,0,floor(t/deltat),omega)
    u4 = cranknicolson(D,deltat,deltax,f,0,0,1000,omega)
    plot2d(u4)
    plot2d(u3)
    //disp(f)
    //disp(u)
    //disp(u2)
    //disp(u3)
endfunction

function a3b(omega,c)
    clf()
    deltax = 0.05
    D = 1
    m = floor(1/deltax)+1
    disp(m)
    f(1) = 0
    f(m) = 0
    for i=2:m-1
        f(i)=sin(%pi*(i-1)*deltax)+c*sin(2*%pi*(i-1)*deltax)
    end
    
    deltat = 0.01
    
        disp(deltat)
    t = 0.07425
    iter = floor(t/deltat)
    u = cranknicolson(D,deltat,deltax,f,0,0,iter,omega)
    plot2d(f)
    plot2d(u)
    t = 0.111375
    u2 = cranknicolson(D,deltat,deltax,f,0,0,floor(t/deltat),omega)
    plot2d(u2)
    t = 0.136125
    u3 = cranknicolson(D,deltat,deltax,f,0,0,floor(t/deltat),omega)
     plot2d(u3)
    //disp(f)
    //disp(u)
    //disp(u2)
    //disp(u3)
endfunction