How to terminate my code?

clc; 
close all; 
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;      

error(1)=1; tol=1e-5;
MI=100; 

u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT); 
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev; 
k=0; 

%while k<MI   %HOW TO USE THIS iterate loop to converge 1 to 0
    for m=1:NT
        for n=1:nodey
            if m>n
                R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n>m
                P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n==m
                Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
            end
        end
    end
    for m=2:NT
        if m==1
            for n=1:nodey
                if n==1   
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
                else
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        else
            for n=1:nodey
                if n==1
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
                else
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        end
        T(:,m)=crank(P,Q,R,S);
        Dt=0.2+Dt;
        t_prev=T;
        
        error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
        if error(m)<tol
            break
        end
    end
    %k=k+1;
   % end
    
   

   
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector


n=length(f);
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
    x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
    fprintf('\n%d = %f\n', i, x(i));
end

[Walter Roberson]

You have the expression

(Dt*/2*DY^2

which has * followed immediately by /

Two lines above that the expression is instead

(Dt/2*DY^2

with no *

[Nathi]

@Walter Roberson My apologise!., I debbugged my code and use updated one.

clc; 
close all; 
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;      

error(1)=1; tol=1e-5;
MI=100; 

u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT); 
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev; 
k=0; 

%while k<MI   %HOW TO USE THIS iterate loop to converge 1 to 0
    for m=1:NT
        for n=1:nodey
            if m>n
                R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n>m
                P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n==m
                Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
            end
        end
    end
    for m=2:NT
        if m==1
            for n=1:nodey
                if n==1   
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
                else
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        else
            for n=1:nodey
                if n==1
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
                else
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        end
        T(:,m)=crank(P,Q,R,S);
        Dt=0.2+Dt;
        t_prev=T;
        
        error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
        if error(m)<tol
            break
        end
    end

Unrecognized function or variable 'crank'.

    %k=k+1;
   % end
    
   

[Walter Roberson]

If you need further assistance, we would need your crank function.

[Nathi]

   
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector


n=length(f);
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
    x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
    fprintf('\n%d = %f\n', i, x(i));
end

[per isakson]

% HOW TO USE THIS iterate loop to converge 1 to 0

What should converge?

[Walter Roberson]

It does converge. But there is a lot of needless output that confuses matters.

clc; 
close all; 
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;      

error(1)=1; tol=1e-5;
MI=100; 

u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT); 
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev; 
k=0; 

%while k<MI   %HOW TO USE THIS iterate loop to converge 1 to 0
    for m=1:NT
        for n=1:nodey
            if m>n
                R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n>m
                P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
            elseif n==m
                Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
            end
        end
    end
    for m=2:NT
        if m==1
            for n=1:nodey
                if n==1   
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
                else
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        else
            for n=1:nodey
                if n==1
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
                else
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        end
        T(:,m)=crank(P,Q,R,S);
        Dt=0.2+Dt;
        t_prev=T;
        
        error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
        error(m)
        if error(m)<tol
            break
        end
    end

1 = 0.000000

2 = 0.000000

3 = 0.000000

4 = 0.000000

5 = 0.000000

6 = 0.000000

7 = 0.000000

8 = 0.000000

9 = 0.000000

10 = 0.000000

11 = 0.000000

12 = 0.000000

13 = 0.000000

14 = 0.000000

15 = 0.000000

16 = 0.000000

17 = 0.000000

18 = 0.000000

19 = 0.000000

20 = 0.000000

21 = 0.000000

22 = 0.000000

23 = 0.000000

24 = 0.000000

25 = 0.000000

26 = 0.000000

27 = 0.000000

28 = 0.000000

29 = 0.000000

30 = 0.000000

31 = 0.000000

32 = 0.000000

33 = 0.000000

34 = 0.000000

35 = 0.000000

36 = 0.000000

37 = 0.000000

38 = 0.000000

39 = 0.000001

40 = 0.000001

41 = 0.000003

42 = 0.000006

43 = 0.000013

44 = 0.000028

45 = 0.000060

46 = 0.000129

47 = 0.000279

48 = 0.000604

49 = 0.001306

50 = 0.002825

51 = 0.006110

52 = 0.013213

53 = 0.028575

54 = 0.061795

55 = 0.133638

56 = 0.289005

ans = 0.2890

1 = 0.000000

2 = 0.000000

3 = 0.000000

4 = 0.000000

5 = 0.000000

6 = 0.000000

7 = 0.000000

8 = 0.000000

9 = 0.000000

10 = 0.000000

11 = 0.000000

12 = 0.000000

13 = 0.000000

14 = 0.000000

15 = 0.000000

16 = 0.000000

17 = 0.000000

18 = 0.000000

19 = 0.000000

20 = 0.000000

21 = 0.000000

22 = 0.000000

23 = 0.000000

24 = 0.000000

25 = 0.000000

26 = 0.000000

27 = 0.000000

28 = 0.000000

29 = 0.000000

30 = 0.000000

31 = 0.000000

32 = 0.000000

33 = 0.000000

34 = 0.000000

35 = 0.000000

36 = 0.000000

37 = 0.000000

38 = 0.000000

39 = 0.000000

40 = 0.000001

41 = 0.000002

42 = 0.000005

43 = 0.000011

44 = 0.000024

45 = 0.000053

46 = 0.000116

47 = 0.000254

48 = 0.000554

49 = 0.001210

50 = 0.002644

51 = 0.005773

52 = 0.012606

53 = 0.027524

54 = 0.060092

55 = 0.131182

56 = 0.286349

ans = 0.0027

1 = 0.000000

2 = 0.000000

3 = 0.000000

4 = 0.000000

5 = 0.000000

6 = 0.000000

7 = 0.000000

8 = 0.000000

9 = 0.000000

10 = 0.000000

11 = 0.000000

12 = 0.000000

13 = 0.000000

14 = 0.000000

15 = 0.000000

16 = 0.000000

17 = 0.000000

18 = 0.000000

19 = 0.000000

20 = 0.000000

21 = 0.000000

22 = 0.000000

23 = 0.000000

24 = 0.000000

25 = 0.000000

26 = 0.000000

27 = 0.000000

28 = 0.000000

29 = 0.000000

30 = 0.000000

31 = 0.000000

32 = 0.000000

33 = 0.000000

34 = 0.000000

35 = 0.000000

36 = 0.000000

37 = 0.000000

38 = 0.000000

39 = 0.000000

40 = 0.000001

41 = 0.000002

42 = 0.000005

43 = 0.000011

44 = 0.000025

45 = 0.000054

46 = 0.000117

47 = 0.000255

48 = 0.000557

49 = 0.001215

50 = 0.002651

51 = 0.005785

52 = 0.012624

53 = 0.027551

54 = 0.060127

55 = 0.131222

56 = 0.286380

ans = 4.0001e-05

1 = 0.000000

2 = 0.000000

3 = 0.000000

4 = 0.000000

5 = 0.000000

6 = 0.000000

7 = 0.000000

8 = 0.000000

9 = 0.000000

10 = 0.000000

11 = 0.000000

12 = 0.000000

13 = 0.000000

14 = 0.000000

15 = 0.000000

16 = 0.000000

17 = 0.000000

18 = 0.000000

19 = 0.000000

20 = 0.000000

21 = 0.000000

22 = 0.000000

23 = 0.000000

24 = 0.000000

25 = 0.000000

26 = 0.000000

27 = 0.000000

28 = 0.000000

29 = 0.000000

30 = 0.000000

31 = 0.000000

32 = 0.000000

33 = 0.000000

34 = 0.000000

35 = 0.000000

36 = 0.000000

37 = 0.000000

38 = 0.000000

39 = 0.000000

40 = 0.000001

41 = 0.000002

42 = 0.000005

43 = 0.000011

44 = 0.000025

45 = 0.000054

46 = 0.000117

47 = 0.000255

48 = 0.000556

49 = 0.001214

50 = 0.002650

51 = 0.005784

52 = 0.012624

53 = 0.027550

54 = 0.060127

55 = 0.131221

56 = 0.286380

ans = 6.9662e-07

    %k=k+1;
   % end
    
      
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector


n=length(f);
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
    x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
    fprintf('\n%d = %f\n', i, x(i));
end
end   

[Nathi]

@Walter RobersonThank you for assit me. Actually, there are two types of convergence here. one is convergence at each time steps and another is convergence at each space steps. In my code, converged at time discretize. That's why i'm using error/tol condition. @per isakson Now I want convergence at each space step i.e., begins with 1 and increase gradually, then diminish gradually to 0 in each 'm'th row. here, my boundary condition are

        if m==1
            for n=1:nodey
                if n==1   %(top boundary condtion i.e (1,1) 1st row and 1st column)
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey  %(bottom boundary condition i.e(nodey,1) last row and 1st col)
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
                else
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        else
            for n=1:nodey
                if n==1    %(at top boundary condition i.e (1,2:NT) 
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey %(at bottom boundary condtion i.e (nodey,2:NT)
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
                else
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        end

If this S get gradually values, then it will store in T. and then if I plot at each column wise in T, it must comes like this "begins with 1 and increase gradually, then diminish gradually to 0 ". For the gradual values and convergence, I'm using "ITERATION".

[Nathi]

@Walter Roberson and @KSSV Please try and include the below condition . It might work.

TT=crank(P,Q,R,S);
T(:,m)=[initial value,TT,endpoint];

t_prev is initial point and t_surf is end point

If get any idea about iteration with boundary condition using while loop., please suggest me.

Answer 1

[KSSV]

You can use for loop instead of while....You know the number of counts right.

clc;
close all;
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;

error(1)=1; tol=1e-5;
MI=100;

u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT);
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev;
E = cell(MI,1) ; 
for  k = 1:MI   %HOW TO USE THIS iterate loop to converge 1 to 0
for m=1:NT
    for n=1:nodey
        if m>n
            R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
        elseif n>m
            P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
        elseif n==m
            Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
        end
    end
end
for m=2:NT
    if m==1
        for n=1:nodey
            if n==1
                S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
            elseif n==nodey
                S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
            else
                S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
            end
        end
    else
        for n=1:nodey
            if n==1
                S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
            elseif n==nodey
                S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
            else
                S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
            end
        end
    end
    T(:,m)=crank(P,Q,R,S);
    Dt=0.2+Dt;
    t_prev=T;

    error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
    if error(m)<tol
        break
    end
    E{k} = error ;
end
%k=k+1;
end

   
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector


n=length(f);
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
    x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
    fprintf('\n%d = %f\n', i, x(i));
end
end

[Nathi]

@KSSV I tried your idea, it's better but not reached with my boundary condition. I'm using iteration for convergence at each space steps and gradual values with my boundary condition.

        if m==1
            for n=1:nodey
                if n==1   %(top boundary condtion i.e (1,1) 1st row and 1st column)
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey  %(bottom boundary condition i.e(nodey,1) last row and 1st col)
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
                else
                    S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        else
            for n=1:nodey
                if n==1    %(at top boundary condition i.e (1,2:NT) 
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
                elseif n==nodey %(at bottom boundary condtion i.e (nodey,2:NT)
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
                else
                    S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
                end
            end
        end

this boundary condtion must satisfy in T. I need the output like If this S get gradually values, then it will store in T. and then if I plot at each column wise in T, it must comes like this "begins with 1 and increase gradually, then diminish gradually to 0 ".