import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def solve_system():
    """
    Solve the system:
        dy/dt - x = 0  --> dy/dt = x
        y + dx/dt = 0  --> dx/dt = -y
    
    This represents a harmonic oscillator: d²y/dt² + y = 0
    """
    
    def derivatives(t, y):
        """
        y[0] = y
        y[1] = x
        
        dy/dt = x
        dx/dt = -y
        """
        dydt = y[1] +     # dy/dt = x
        dxdt = -y[0]     # dx/dt = -y
        return [dydt, dxdt]
    
    # Initial conditions: y(0) = 1, x(0) = 0
    y0 = [1.0, 0.0]  # [y(0), x(0)]
    
    # Time span
    t_span = (0, 10)
    t_eval = np.linspace(0, 10, 200)
    
    # Solve the system
    sol = solve_ivp(
        fun=derivatives,
        t_span=t_span,
        y0=y0,
        t_eval=t_eval,
        method='RK45'  # 4th-order Runge-Kutta
    )
    
    return sol.t, sol.y[0], sol.y[1]

# Solve
t, y, x = solve_system()

# Plot results
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(t, y, 'b-', label='y(t)')
plt.plot(t, x, 'r--', label='x(t)')
plt.xlabel('t')
plt.ylabel('y, x')
plt.title('Solution: dy/dt = x, dx/dt = -y')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(y, x, 'g-')
plt.xlabel('y')
plt.ylabel('x')
plt.title('Phase Plot (y vs x)')
plt.grid(True)
plt.axis('equal')

plt.tight_layout()
plt.savefig('solution_plot.png', dpi=150)
plt.show()

# Print some values
print("Sample results:")
print(f"{'t':>8} {'y(t)':>12} {'x(t)':>12}")
print("-" * 36)
for i in range(0, len(t), len(t)//5):
    print(f"{t[i]:>8.2f} {y[i]:>12.6f} {x[i]:>12.6f}")