運動方程式を解きたいです。強制外力項に変位が入って計算できません。

Question

import sympy as sp
import math
import cv2
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import glob
import pandas as pd

options = {}
options['strings_to_formulas'] = False
options['strings_to_urls'] = False
with pd.ExcelWriter('eq_CL=tan.xlsx', mode='w') as writer:
    results = []
    for i in range(48) :
        print(i)
        Vr = i*0.5 + 0.5
        print(Vr)
        π = math.pi
        m = 0.0156
        δ = 0.0356167639226128
        D = 0.016
        H = 0.04
        L = 0.3
        ρ = 1.18
        M = (m/(ρ*H*D*L))*(D/H)
        
        lists = []
        sp.init_printing()
        sp.var('t')
        y = sp.Function('y')(t) #xはtの関数(Function)
        α = (4*π*δ)/Vr
        ωn = (2*π)/Vr
        a = ((-1)/(15*(180/π)))
        CL = a*sp.atan(y/Vr)
        F = CL/(2*M)
        T = Vr #=2π/ωn
        eq2 = sp.Eq( sp.diff(y, t ,2)+α*sp.diff(y, t, 1)+(ωn*ωn)*y,  F )　#運動方程式
        yre = sp.dsolve(eq2, ics={y.subs(t,0):0, y.diff(t,1).subs(t,(T/4)):0})
        
        for j in range(0, 4001):
            print(j)
            lists.append(yre.subs(t, j*0.0005))
        results.append(lists)

    df = pd.DataFrame(results).T
    
    df.to_excel(writer)

上記のようなコードで，ばねマスダンパ系の強制振動を解きたいです．
強制外力がFとしており，FをCLとMで構成しています．
CLに定数を与えるとうまく実行できます．

今回したいことですが、
CLに求める変位yを入れることです．
もちろん，求めるｙを使用してｙを求めようとしているため計算ができません．（当たり前なのですが）

エクセルのゴーシーク機能のように解を収束させる形でこの方程式を解くことはできないでしょうか？

よろしくお願いいたします．