Vectorization for temperature simulation

I'm new to Matlab and I would like to know tips in order to reduce the processing time of such program. The scheme of the actual code is very similar to typical C++ code and I would like to vectorize the for and if loops in order to decrease the processing time. The code is added just to make you know the overall idea of what I would like to obtain.

close all             %Tanca totes les finestres obertes
clear all             %Removes variables from memory
clc                   %Clears command window
tic                   %inicia comptador temps simulació

trans=fopen('trans_matlab.txt','wt') 

%Physical properties and geometrical data input
Tair=25;
Tinitial=100;
h=10;
k=5;
qv=500;
thick1=0.05;
thick2=0.1;
density=1000;
Cp=800;

%Parameters imput
max_error=0.000001;   %Error màxim a partir del cual no hi ha canvis en els resultats
max_iter=1000000;
r=0.9;
N=10;                 %Nodes paret1          
M=10;                 %Nodes paret2
Inc_x1=thick1/N;
Inc_x2=thick2/M;
num_t=23000;          %temps de realització
Inc_t=20;

%Initial map temperatures
t=1;
for i=1:(N+M+1)       %Hipòtesi inicial: tots els nodes tenen T i Ta = Tair                 
    T(i,t)=Tinitial;      %(N+M+1) pq no deixa començar a 0  
    Ta(i,t)=Tinitial;
end

%GAUS SEIDEL
while t<num_t
    %New instant
    t0=t;
    t=t+1;    
    for i=1:(N+M+1)
        Ta(i,t)=T(i,t0);
    end
    %Temperature calculation
    for j=0:max_iter                  
        for i=1:(N+M+1)
            if i==1
            b=0;
            c=k/Inc_x1;
            d=qv*Inc_x1/2;
            e=density*Cp*Inc_x1/(2*Inc_t);
            a=b+c+e;            
            T(i,t)=(c.*Ta(i+1,t)+d+e.*T(i,t0))/a;
            end
            if (i>1)&(i<(N+1))        
            b=k/Inc_x1;
            c=k/Inc_x1;
            d=qv*Inc_x1;
            e=density*Cp*Inc_x1/(Inc_t);            
            a=b+c+e;
            T(i,t)=(b.*Ta(i-1,t)+c.*Ta(i+1,t)+d+e.*T(i,t0))/a;      
            end
            if i==(N+1)
            b=k/Inc_x1;
            c=k/Inc_x2;
            d=qv*Inc_x1/2;
            e=density*Cp*Inc_x1/(2*Inc_t)+density*Cp*Inc_x2/(2*Inc_t);
            a=b+c+e;
            T(i,t)=(b.*Ta(i-1,t)+c.*Ta(i+1,t)+d+e.*T(i,t0))/a;
            end
            if (i>(N+1))&(i<(N+M+1))
            b=k/Inc_x2;
            c=k/Inc_x2;
            e=density*Cp*Inc_x2/(Inc_t);
            a=b+c+e;
            T(i,t)=(b.*Ta(i-1,t)+c.*Ta(i+1,t)+e.*T(i,t0))/a;
            end
            if i==(N+M+1)
            b=k/Inc_x2;
            c=h;
            e=density*Cp*Inc_x2/(Inc_t*2);
            a=b+c+e;
            T(i,t)=(b.*Ta(i-1,t)+h*Tair+e.*T(i,t0))/a;
            end
            %Error calculation-trobar l'error maxim de la matriu error(i,t)
            errort=0;
            error(i,t)=abs(T(i,t)-Ta(i,t));
            if error(i,t)>errort
            errort=error(i,t);
            else
            errort=errort;
            end
        end
        if errort>max_error
            for i=1:(N+M+1)
                Ta(i,t)=T(i,t).*r+Ta(i,t).*(1-r);
            end
        else
            break
        end
    end
    iter(t)=j;
end

toc
%Print Results
t=1;
fprintf(trans,'time\t T1\t TN+1\t TN+M+1\t iter\n');

while(t<=num_t)
time=t*Inc_t;
fprintf(trans,'%.2f\t %.2f\t %.2f\t %.2f\t %.2f\n',time,T(1,t),T(N+1,t),T(N+N+1,t),iter(t));
t=t+1;
end
fclose(trans);