# -*- coding: utf-8 -*-
"""
Created on Sun Oct  9 11:00:22 2022

@author: Chris
Lab 10 Q2
"""

import numpy as np
import matplotlib.pyplot as plt

def f(x,y):
    f = x*(1+4*y**2)
    return f


def rk2solve(f,x0,y0,xmax,h):
    xvals = np.array([x0])
    yvals = np.array([y0])
    
    x = x0
    y = y0
    while x<xmax:
        
        k1 = h*f(x,y)
        k2 = h*f(x+0.5*h,y+0.5*k1)
        yn = y + k2
        
        xn = x + h
    
        xvals = np.append(xvals,xn)
        yvals = np.append(yvals,yn)
        x,y = xn,yn
        #print(xn,yn)
        
    return xvals,yvals
    
    
x0 = 0
y0 = 0
h = 0.01

xmax = 1.0

xvals,yvals = rk2solve(f,x0,y0,xmax,h)
plt.plot(xvals,yvals)