Would interp1 be accurate enough for you? If so, try this:
a = 5;
n = 1:4;
theta_teta = -10;
delta_theta = 40;
theta = -10:1:30;
for j = 1: length (n);
for i = 1:length(theta)
f(i) = 1-exp(-a*((theta(i)-theta_teta)/delta_theta)^n(j));
end
plot (theta, f, 'LineWidth', 2);
max_f = sum(f);
fprintf('For n = %d, the sum of f (for some reason called max_f) = %f.\n', j, max_f);
grid on;
hold on;
% Find x value for y = 0.5
x50(j) = interp1(f, theta, 0.5);
line([x50(j), x50(j)], [0, 0.5], 'LineWidth', 2, 'Color', 'k');
fprintf('For n = %d, y = 0.5 at x = %f.\n', j, x50(j));
% Find x value for y = 0.9
x90(j) = interp1(f, theta, 0.9);
line([x90(j), x90(j)], [0, 0.9], 'LineWidth', 2, 'Color', 'k');
fprintf('For n = %d, y = 0.9 at x = %f.\n', j, x90(j));
end
legend('n=1', 'n=2', 'n=3', 'n=4', 'Location', 'east');
% Put lines across at y = 0.9 and 0.5
line(xlim, [0.5, 0.5], 'LineWidth', 2, 'Color', 'k');
line(xlim, [0.9, 0.9], 'LineWidth', 2, 'Color', 'k');
% Label the axes
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);
For n = 1, the sum of f (for some reason called max_f) = 32.540191.
For n = 1, y = 0.5 at x = -4.439359.
For n = 1, y = 0.9 at x = 8.435958.
For n = 2, the sum of f (for some reason called max_f) = 24.668279.
For n = 2, y = 0.5 at x = 4.894287.
For n = 2, y = 0.9 at x = 17.153014.
For n = 3, the sum of f (for some reason called max_f) = 19.624252.
For n = 3, y = 0.5 at x = 10.702175.
For n = 3, y = 0.9 at x = 20.896785.
For n = 4, the sum of f (for some reason called max_f) = 16.262926.
For n = 4, y = 0.5 at x = 14.406551.
For n = 4, y = 0.9 at x = 22.955276.
