function pumps(v1,v2,n,t)
    R = 2 + (5-2) * rand(1,n);
    T = -log(rand(1,n))/0.5;
    S1 = zeros(1,n);S2 = zeros(1,n);W1 = zeros(1,n);W2 = zeros(1,n);W1n = zeros(1,n);W2n = zeros(1,n);
    S1(1) = R(1)/v1; S2(1) = R(1)/v2; 
    for i = 2:n
        S1(i) = R(i)/v1;
        S2(i) = R(i)/v2;
        W1(i) = max(W1(i-1) + S1(i-1) - T(i),0);
        W2(i) = max(W2(i-1) + S2(i-1)-T(i),0);        
        for q = 1:i
            W1n(i) = W1(q) + W1n(i);
            W2n(i) = W2(q) + W2n(i);
        end%for
        W1n(i) = W1n(i) / i;
        W2n(i) = W2n(i) / i;
    end%for
    figure
    plot(1:n,W1n,1:n,W2n);
    hleg = legend('V=1','V=5');
    set(hleg,'Location','NorthWest');
    top = title('Average Wait time for varying pump capacties');
    y = ylabel('Average Wait Time'); x = xlabel('Customer');
    
  %Work in System  
    
  arrival = zeros(1,n);aw1 = zeros(1,n);aw2 = zeros(1,n);
  arrival(1) = T(1); aw(1) = arrival(1);
  for x =2:n
      arrival(x) = arrival(x-1) + T(x);
      aw1(x) = arrival(x) + W1(x);
      aw2(x) = arrival(x) + W2(x);
  end
  L1 = zeros(1,t); L2 = zeros(1,t);
  for j = 1:t
    u = 1;
    while u<n+1 && arrival(u) < j
      if aw1(u) > j         
          L1(j) = S1(u) + L1(j);         
      else%if
        if(aw1(u)+S1(u) > j)
          timediff = j - aw1(u);
          if timediff*v1 < S1(u)          
            workdone = max(S1(u) - timediff*v1,0);
            L1(j) = L1(j) + workdone;            
          end%if
        end%if
      end%if
      u = u + 1;      
    end%while 
  end%for
 
    figure
    plot(1:t,L1); 
    hleg = legend('V=1');
    set(hleg,'Location','NorthWest');
    top = title('Work to be done for V=1');
    y = ylabel('Work to be done'); x = xlabel('Customer');
    
  
  
   for j = 1:t
    u = 1;
    while u<n+1 && arrival(u) < j
      if aw2(u) > j         
          L2(j) = S2(u) + L2(j);         
      else%if
        if(aw2(u)+S2(u) > j)
          timediff = j - aw2(u);
          if timediff*v1 < S2(u)          
            workdone = max(S2(u) - timediff*v2,0);
            L2(j) = L2(j) + workdone;            
          end%if
        end%if
      end%if
      u = u + 1;      
    end%while 
  end%for
  
    figure
    plot(1:t,L2); 
    hleg = legend('V=5');
    set(hleg,'Location','NorthWest');
    top = title('Work to be done for V=5');
    y = ylabel('Work to be done'); x = xlabel('Customer');
    
  
  
end%function