How to program difference equations over time?

Hey,

can somebody refer me to a good tutorial on how to program difference equations over time? I'm currently trying to program something like this, but I'm stuck with step 3):

1. Set initial values for some variables which I want to track over time (here: uI, uM, Omega and Equ)
2. Define some other parameters and helpful equations
3. Solve difference equations for 3 other variables
4. The variables from 3) + equations for the law of motion of the variables from 1), define the new starting values
5. Start over with 1)

I only found info on differential equations and played around a bit with tspan and a (t,x) function, but unfortunately it didn't yield anything useful. So I'm looking for a suitable function to solve step %Calculate Lambda and Theta% in the code and make it dynamic over time:

if true
  % % Variables for Difference equations
%
thet_I = x(1,1);                                                            
thet_M = x(2,1);                                                            
lambd = x(3,1);                                                             
%
% Initial values for u, Omega and Eq
%
uM_0 = 0.0703;
uI_0 = 0.0689;
Omeg_0 = 0.147;
Equ_0 = 0.0068;
%
% set values for Parameters of the model
% 
alph = 0.36;                                              
bet_B = 0.98895;                                          
bet_H = 0.9956;                                           
delt = 0.005354932;                                       
gamm = 0.1121;                                            
z_I = 1.0704;                                             
z_M = 1;                                                  
b = 0.4;                                                  
et = 0.5;                                                 
kapp = 0.213;                                             
Lbar_I = 0.66;                                            
Lbar_M = 1;                                               
s_M = 0.03451;                                            
s_I = 0.03451;                                            
mu = 0.4430;                                              
iot = 0.72;
T = 100;                                                  %Number of time periods%
%
% Define variables for Difference equation
%
rh_H = (1/bet_H)-1;                                                           
rh_B = (1/bet_B)-1;                                            
lambd = (delt + rh_B + 1)/gamm;                                
r_M = delt + rh_H;                                             
r_I = r_M + (gamm * (rh_B-rh_H))/(1+delt+rh_B);                
bet_FM = 1/(1+r_M);                                            
bet_FI = 1/(1+r_I);                                            
q_M = mu*thet_M^(-iot);                                        
q_I = mu*thet_I^(-iot);                                        
pi_M = thet_M*q_M;                                             
pi_I = thet_I*q_I;                                             
L_M = (1-uM(time,1)) * Lbar_M;                                 
L_I = (1-uI(time,1)) * Lbar_I;                                 
K_I = ((alph*z_I)/r_I)^(alph-1)*L_I;                           
K_M = ((alph*z_M)/r_M)^(alph-1)*L_M;                           
Outp_M = (1/L_M)*(z_M * K_M^alph * L_M^(1-alph));              
Outp_I = (1/L_I)*(z_I * K_I^alph * L_I^(1-alph));              
w_M=(et*(Outp_M-r_M*(K_M/L_M)+kapp*thet_M))/(1-(1-et)*b);     
w_I=(et*(Outp_I-r_I*(K_I/L_I)+kapp*thet_I))/(1-(1-et)*b);      
KT = K_I + K_M;                                                
R_B = gamm*lambd-delt;
%
% Calculate level of Lambda and Theta
%
eqs = ones(length(x),1);                                       
eqs(1,1) = (kapp/q_I)*(1-bet_FI*(1-s_I))-bet_FI*(Outp_I-r_I*(K_I/L_I)-w_I);  
eqs(2,1) = (kapp/q_M)*(1-bet_FM*(1-s_M))-bet_FM*(Outp_M-r_M*(K_M/L_M)-w_M); 
eqs(3,1) = (1-(1/lambd))*alph*z_M*((Omeg(time+1,1)+Equ(time+1,1)-lambd*Equ(time+1,1))/L_M)^(alph-1)-(1/lambd)+gamm-alph*z_I*(lambd*Equ(time+1,1)/L_I)^(alph-1);
%
%  Initialize time series 
%  
uI = zeros(T,1);
uM = zeros(T,1);
Omeg = zeros(T,1);
Equ = zeros(T,1);
uM(1) = uM_0;
uI(1) = uI_0;
Omeg(1) = Omeg_0;
Equ(1) = Equ_0;
%
%  Simulate model
%
for time = 1:T-1
  uI(time+1,1) = uI(time,1)+s_I*(1-uI(time,1))-pi_I*uI(time,1);
  uM(time+1,1) = uM(time,1)+s_M*(1-uM(time,1))-pi_I*uM(time,1);
  Omeg(time+1,1) = bet_H*(1+r_M-delt)*Omeg(time,1);
  Equ(time+1,1) = bet_B*R_B*Equ(time,1);
end;
%
%  plot graphs
%
figure(1)
plot((1:(1+T-1)),uI)
xlabel('year');
ylabel('Unemployment rate I');
title('Development of unemployment');
end