% Define the parameters
lambda = 8;
mu = 10;
alpha = 50;
beta = 20;
c = 1;
h = 0.05;

% Define the expected cost function
function cost = expected_cost(n)
  rho = lambda / mu;
  pi0 = 1 / (sum(rho .^ (0:n) ./ factorial(0:n)) + rho .^ (n+1) / (factorial(n) * (1 - rho)));
  pi = rho .^ (0:n) .* pi0 / factorial(0:n);
  Lq = sum((1:n) .* pi(2:end));
  Wq = Lq / (mu - lambda);
  W = Wq + 1 / mu;
  cost = alpha + beta * w + c * h * (n - Lq);
endfunction

% Define the range of values for n
n_range = 1:50;

% Compute the expected cost for each value of n
costs = arrayfun(@expected_cost, n_range);

% Plot the results in 3D
mesh(n_range, ones(size(n_range)), costs);
xlabel('Storage capacity, n');
ylabel('Dummy variable');
zlabel('Expected cost');

