Matlab code for solve a differential equation with Alternating Direction Implicit (ADI)

      

        %ADI Summary of this function goes here
        %  Detailed explanation goes here
      
 n=98;
        ll=0;lr=0;lt=1;lb=1;
        mallatotal=zeros(n+2,n+2);
        mallatotal(:,1)=ll;
        mallatotal(:,n+2)=lr;
        mallatotal(1,:)=(lb * (4/n)*(1/(n+2))* ((n/2+.5)^2-( (1:n+2)-(n+3)/2).^2 ))

        mallatotal(n+2,:)=(lt * (4/n)*(1/(n+2))* ((n/2+.5)^2-( (1:n+2)-(n+3)/2).^2 ))

        %esquinas

        mallatotal(1,1)=0*(lb+ll)/2;

        mallatotal(1,n+2)=0*(lb+lr)/2;

        mallatotal(n+2,1)=0*(lt+ll)/2;

        mallatotal(n+2,n+2)=0*(lt+lr)/2;

        rh=zeros(n,1);

        h=2/(n+1);

        k=h/2;r=k/h^2;ro=2/r;

        vold=zeros(n);vnew=vold;vint=vold;

        sub=-ones(n-1,1);

        d=(ro+2)*ones(n,1);

        id2=diag(sub,-1)+diag(d)+diag(sub,1);

        sub=-ones(n-1,1);

        d=(ro+2)*ones(n,1);

        id1=diag(sub,-1)+diag(d)+diag(sub,1);

        normaold=sqrt(sum(diag(vnew'*vnew)));

        normanew=10000000+.5*normaold; %cualquier valor que las haga suficientemente diferentes

        t=0;

        while abs(normanew-normaold)/normanew>.00001

        for j=1:n

        if j==1

        rh(:)=mallatotal(2:n+1,1)+(ro-2)*vold(:,1)+vold(:,2);

        elseif j==n

        rh(:)=vold(:,n-1)+(ro-2)*vold(:,n)+mallatotal(2:n+1,n+2);

        else

        rh(:)=vold(:,j-1)+(ro-2)*vold(:,j)+vold(:,j+1);

        end

        rh(1)=rh(1)+1*mallatotal(1,1+j);rh(n)=rh(n)+1*mallatotal(n+2,1+j);

        v=id1\rh;

        vint(:,j)=v(:);

        end

        for i=1:n

        if i==1

        rh(:)=mallatotal(1,2:n+1)+(ro-2)*vint(1,:)+vint(2,:);

        elseif i==n

        rh(:)=vint(n-1,:)+(ro-2)*vint(n,:)+mallatotal(n+2,2:n+1);

        else

        rh(:)=vint(i-1,:)+(ro-2)*vint(i,:)+vint(i+1,:);

        end

        rh(1)+1*mallatotal(1+j,1);rh(n)=rh(n)+1*mallatotal(n+2,1+j);

        v=id2\rh;

        vnew(i,:)=v(:)';

        end

        t=t+1;

        vold=vnew;

        normaold=normanew;

        normanew=sqrt(sum(diag(vnew'*vnew)));

        mesh(vnew);

        pause(.1)

        end

        t

        normanew;

        vnew;

        mallatotal(2:n+1,2:n+1)=vnew(1:n,1:n);

        mesh(mallatotal);

        AXIS([0 n+3 0 n+3 0 1]);