clc
clear
close all

%% Input and output data
V = (-100:5:-65)';
Iss = [-10; -8; -7.5; -6.5; -4; -2; 0.5; 7];

%% Scale voltage for numerical stability
V_mean = mean(V);
V_std = std(V);
z = (V - V_mean) / V_std;

%% Design matrices
X1 = [ones(size(z)), z];
X2 = [ones(size(z)), z, z.^2];
X3 = [ones(size(z)), z, z.^2, z.^3];

%% Least-squares parameter estimation
theta1 = X1 \ Iss;
theta2 = X2 \ Iss;
theta3 = X3 \ Iss;

%% Model outputs at the measured data points
I1_hat = X1 * theta1;
I2_hat = X2 * theta2;
I3_hat = X3 * theta3;

%% Mean squared errors
MSE1 = mean((Iss - I1_hat).^2);
MSE2 = mean((Iss - I2_hat).^2);
MSE3 = mean((Iss - I3_hat).^2);

%% Display parameters and errors
fprintf('Voltage scaling: z = (V - %.6f) / %.6f\n\n', V_mean, V_std)

fprintf('First-degree model parameters [theta0; theta1]:\n')
disp(theta1)
fprintf('First-degree MSE = %.6f\n\n', MSE1)

fprintf('Second-degree model parameters [theta0; theta1; theta2]:\n')
disp(theta2)
fprintf('Second-degree MSE = %.6f\n\n', MSE2)

fprintf('Third-degree model parameters [theta0; theta1; theta2; theta3]:\n')
disp(theta3)
fprintf('Third-degree MSE = %.6f\n\n', MSE3)

%% Smooth curves for plotting
V_plot = linspace(min(V), max(V), 500)';
z_plot = (V_plot - V_mean) / V_std;

X1_plot = [ones(size(z_plot)), z_plot];
X2_plot = [ones(size(z_plot)), z_plot, z_plot.^2];
X3_plot = [ones(size(z_plot)), z_plot, z_plot.^2, z_plot.^3];

I1_plot = X1_plot * theta1;
I2_plot = X2_plot * theta2;
I3_plot = X3_plot * theta3;

%% Four subplots in one figure
figure

subplot(2,2,1)
plot(V, Iss, 'ko', 'MarkerFaceColor', 'k')
hold on
plot(V_plot, I1_plot, 'b-', 'LineWidth', 1.5)
grid on
xlabel('Voltage (V)')
ylabel('Current (I)')
title(sprintf('First-Degree Model, MSE = %.4f', MSE1))
legend('Measured Data', 'Fitted Model', 'Location', 'northeast')

subplot(2,2,2)
plot(V, Iss, 'ko', 'MarkerFaceColor', 'k')
hold on
plot(V_plot, I2_plot, 'r-', 'LineWidth', 1.5)
grid on
xlabel('Voltage (V)')
ylabel('Current (I)')
title(sprintf('Second-Degree Model, MSE = %.4f', MSE2))
legend('Measured Data', 'Fitted Model', 'Location', 'northeast')

subplot(2,2,3)
plot(V, Iss, 'ko', 'MarkerFaceColor', 'k')
hold on
plot(V_plot, I3_plot, 'g-', 'LineWidth', 1.5)
grid on
xlabel('Voltage (V)')
ylabel('Current (I)')
title(sprintf('Third-Degree Model, MSE = %.4f', MSE3))
legend('Measured Data', 'Fitted Model', 'Location', 'northeast')

subplot(2,2,4)
plot(V, Iss, 'ko', 'MarkerFaceColor', 'k')
hold on
plot(V_plot, I1_plot, 'b-', 'LineWidth', 1.5)
plot(V_plot, I2_plot, 'r-', 'LineWidth', 1.5)
plot(V_plot, I3_plot, 'g-', 'LineWidth', 1.5)
grid on
xlabel('Voltage (V)')
ylabel('Current (I)')
title('Comparison of Polynomial Models')
legend('Measured Data', 'Degree 1', 'Degree 2', 'Degree 3', ...
       'Location', 'northeast')