# -*- 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

x0 = 0
y0 = 0
h = 0.01

xmax = 1.0

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)
    
plt.plot(xvals,yvals)