Dos gráficas en una

buenos días a todos soy Gustavo estoy dando mis primeros pasos en python, bueno tengo algo de experiencia con lo básico.
El problema que tengo y espero que algunos de ustedes le agradecería que me pueda ayuda es el siguiente:
Estoy intentando crear un bucle donde me permita visualizar la gráfica de dos soluciones de una EDO lineal, la numérica y la exacta.
El código fuente que estoy tratando de ejecutar es el siguiente:
“”" ********************************************************
*** Código fuente donde se aplica el método de Euler ***
*** hacia adelante para resolver una EDO lineal de ***
*** primer orden. ***
******************************************************** “”"

librería

import os # permite limpiar pantalla.
import numpy as np
import matplotlib.pyplot as plt
from scipy import special # calcula la función error.

“”" ************************************
*** Definición de las variables. ***
************************************
a → Condición inicial para la solución.
x0,y0 → Condiciones iniciales en x e y.
h → Separación de los intervalos.
xx → Valor de x para conocer y. “”"

os.system (“cls”) # limpia pantalla
x0 = float(input("\n Introducir el valor de "x": “))
a = float(input(”\n Introducir el valor de "y": “))
xx = float(input(”\n Introducir el valor de "x" requerido para conocer "y": “))
h = float(input(”\n Introducir la separación de los intervalos: "))

“”" ************************************
*** Definición de las funciones. ***
************************************ “”"
def SolucionExacta(x):
# Solución de la EDO lineal.
return (4.0 * a * np.exp(x2) + 2 * x + np.sqrt(np.pi) * np.exp(x2) *
special.erf(x))/4

def funcion(t,y):
return 2 * t * y - t**2 + 1.0 # dy/dt EDO lineal

Realización de los cálculos

y0 = a # condiciones iniciales
y = y0
x = x0
xh = 0
iter = 0
while xh <= xx:
X = np.arange(0,xh,0.005)
y += h * funcion(xh,y) # solución numérica de la EDO lineal.
exacta = SolucionExacta(xh) # solución exacta de la EDO lineal.
plt.figure()
plt.plot(X, y, color = ‘r’, label = ‘Sol.num.’)
plt.plot(X, exacta, color = ‘g’, label = ‘Sol.exa.’)
xh += h
plt.xlabel(r’$x$‘)
plt.xlabel(r’$y(x)$')
plt.show()


Dentro del cuerpo del bucle while, es donde pretendo que las dos gráficas se desplieguen (no es una animación).
Al ejecutar la aplicación, pueden colocar como parámetros de entradas:
x0 = 0
a = 1.5
xx = 3
h = 0.10

Le agradecería de antemano cualquier ayuda que puedan brindarme para superar este inconveniente.

No estoy seguro, pero ahi va:

import os  # permite limpiar pantalla.
import numpy as np
import matplotlib.pyplot as plt
from scipy import special  # calcula la función error.

""" ************************************
*** Definición de las variables. ***
************************************
a → Condición inicial para la solución.
x0,y0 → Condiciones iniciales en x e y.
h → Separación de los intervalos.
xx → Valor de x para conocer y. """

os.system("clear")  # limpia pantalla
x0 = float(input('\n Introducir el valor de "x": '))
a = float(input('\n Introducir el valor de "y": '))
xx = float(input('\n Introducir el valor de "x" requerido para conocer "y": '))
h = float(input('\n Introducir la separación de los intervalos: '))


""" ************************************
*** Definición de las funciones. ***
************************************ """


def SolucionExacta(x):
    # Solución de la EDO lineal.
    return (4.0 * a * np.exp(x) + 2 * x + np.sqrt(np.pi) * np.exp(x) * special.erf(x))/4


def funcion(t, y):
    return 2 * t * y - t**2 + 1.0  # dy/dt EDO lineal


y0 = a  # condiciones iniciales
y = y0
x = x0
xh = 0
iter = 0

ys = []
exactas = []

while xh < xx:
    ys.append(h * funcion(xh, y))
    exactas.append(SolucionExacta(xh))
    xh += h

X = np.arange(0, xx, h)
plt.figure()
plt.plot(X, ys, color='r', label='Sol.num.')
plt.plot(X, exactas, color='g', label='Sol.exa.')
plt.xlabel(r'$x$')
plt.xlabel(r'$y(x)$')
plt.show()

Figure_1