% Plot the isosurfaces of I1 and I2

function [] = isosurfaces(y0, I)

    I1 = @(u) norm(u)^2;
    I2 = @(u) (u(1)^2/I(1) + u(2)^2/I(2) + u(3)^2/I(3))/2;

    Np = 50;
    x = linspace(-1.5, 1.5, Np);
    y = linspace(-1.5, 1.5, Np);
    z = linspace(-1.5, 1.5, Np);
    [X,Y,Z] = meshgrid(x, y, z);

    I1val = zeros(size(X));
    I2val = zeros(size(X));
    for i = 1:Np
        for j = 1:Np
            for k = 1:Np
                I1val(i, j, k) = I1([X(i, j, k), Y(i, j, k), Z(i, j, k)]);
                I2val(i, j, k) = I2([X(i, j, k), Y(i, j, k), Z(i, j, k)]);
            end
        end
    end

    isovalue1 = I1(y0);
    surf1 = isosurface(X, Y, Z, I1val, isovalue1);
    p1 = patch(surf1);
    isonormals(x, y, z, I1val, p1);
    set(p1, 'FaceColor', 'red', 'EdgeColor', 'none', 'FaceAlpha', 0.1); % set color, mesh and transparency level of the surface
    
    daspect([1, 1, 1]);
    view(3); 
    camlight;
    lighting gouraud;

    isovalue2 = I2(y0);
    surf2 = isosurface(X, Y, Z, I2val, isovalue2);
    p2 = patch(surf2);
    isonormals(x, y, z, I2val, p2);
    set(p2, 'FaceColor', 'blue', 'EdgeColor', 'none', 'FaceAlpha', 0.1); % set color, mesh and transparency level of the surface

end