Analyzing a constant current source
Quick analysis of this circuit

First, one must do some assumptions which stems from the art of electronic engineering.
One "verify" afterward that there no contradiction between the results and the assumptions.  

These assumptions are :   


The current in the base of the transistor is so low that the voltage at the point common between rp1 and rp2 is the voltage without load.   
Vbe = 0.7 Volts  
Iemitter = Icollector


Then :  


Vbase = VBat ( rp1 + rp2)/(rp1 + rp2+ rs) = 3.6 2.7 Volts (the Zener has no effect, see my comment)    
Vemitter = 3.6 2.7 - 0.7 = 2.9 2.0 Volts  
Iemitter = 2.9 2.0/Re =  5.8 4 mA  
Icollector = 5.8 4 mA  


Let's see if there are no contradiction with the initial assumptions :  


Ibase = ICollector/Beta =
   Icollector /(Something more than 40 according to the datasheet of the transistor)
   &lt; 0.145 0.100 mA
Ibase comes from some equivalent resistor which is less than 750 Ohms. This implies that the voltage drop is no more than 750 Omhs * 0.145 0.100 mA  = 0.108 0.075 Volts
0.108 0.075 Volts can be considered as negligible compared with 0.7 Volts.    
Concerning Vbe. This voltage is the voltage that would develop a diode. This voltage is  0.12 Log[10,Idiode / 10^-10] = 0.7 Volts. Note that this value is very insensitive to Idiode because of the Log  
Iemitter = Icollector is true if Vce > 1 Volt (ie the transistor works in is active region). Here we have Vce = VBatt - 100 5.8 4 mA -  2.9 2.0 = 5.5 6.6 Volts    


Solution of the problem with the Wolfram Language (Mathematica)  

rp1=250;
rp2=750;
rs=1500;
rl=100;
VBat=9;
re=500;
beta=40;  

eqns = {
   VBat - ic rl - vce - ie re == 0 ,
   VBat - i3 rs - i3 rp1 - i5 rp2 == 0,
   ie == ib + ic,
   i1 == ic + i3,
  (* i1 == i5 + ie, redundant with the line just above *)
   i3 == i5 + ib,
   beta == ic / ib,
   i5 rp2 == ie re + 0.7 (* new *)
   };

res00=Solve[eqns,{ic,vce,ie,i3,i5,i1,ib}][[1]];
res00  //Column  




solution with a more refined model of the jonction base-emitter (Vbe~0.7Volts)  

eqns01 = {
   VBat - ic rl - vce - ie re == 0 ,
   VBat - i3 rs - i3 rp1 - i5 rp2 == 0,
   ie == ib + ic,
   i1 == ic + i3,
  (* i1 == i5 + ie, redundant with the line just above *)
   i3 == i5 + ib,
   beta == ic / ib,
   i5 rp2 == ie re + 0.12 Log[10,ib / 10^-10] (* a more refined model  *)
   };


res01=FindRoot[Evaluate @ eqns01,Evaluate[List @@@ res00]]
res01 //Column