restart; diffeq := collect({diff(y(z),z,z)-(-2*z^5+4*z^3+z^4*a-2*z-a)*(diff(y(z),z))/((z-1)^3*(z+1)^3)-(-z^2*b+(-c-2*a)*z-d)*y(z)/((z-1)^3*(z+1)^3),y(0)=1,D(y)(0)=0},diff,factor);

PCUsKC1JJWRpZmZHJSpwcm90ZWN0ZWRHNiQtRiU2JC1JInlHNiI2I0kiekdGLEYuRi4iIiIqKiwqKiRGLiIiJCEiIyomRi4iIiNJImFHRixGL0YvRi5GNkY3Ri9GLywmRi5GL0YvRi9GNCwmRi5GLyEiIkYvRjRGKEYvRjoqKiwqKiZGLkY2SSJiR0YsRi9GLyomRi5GL0kiY0dGLEYvRi8qJkYuRi9GN0YvRjZJImRHRixGL0YvRipGL0Y5ISIkRjhGQ0YvLy1GKzYjIiIhRi8vLS1JIkRHNiRGJkkoX3N5c2xpYkdGLDYjRitGRkZH

a, b, c, d := 1, 1/3, 1/2, 3;

NiYiIiIjRiMiIiQjRiMiIiNGJQ==

evalf[51](HeunD(a, b, c, d, 1/3));

JCJUcD8wLnVXJlIrNSNlU0p5Kz1SRCZSY1p1OlA3ISNd

myHeunD := diffeqtoproc(diffeq, y(z)): myHeunD(1/3, 50);

JCJUdj8wLnVXJlIrNSNlU0p5Kz1SRCZSY1p1OlA3ISNd

myHeunD(1/3, 2000);

JCJcaHI0MyFvXUJ4IzNEOylHV29xKioqPk9MZzZmZysxZFpkWVpMP0xMJ2V2USFwOzVlTDkjR0d1Tnhud0JqXGtIdGVoIjRCbGlDMzNRXyw8J1wuQSFvWkUzXz9uOTt5J1IzLyZwJHp1LHJqI1JmOFp0aWFCQlAuLSx1Jm82ZycpUiVbUzojcD0mPiwoR0ZTKDNzIjRHLHElKVEpeSpbKlxxRyoqZjtHXXIvJEh0VGFDK2dBXjNMODMzLkxJT3NBJVxdeEQnKSlbcyZIKFtULigqPjxEPHBCdnhCdXVBLChmJ3BDOzUtRzcsczJQNXYkKkc0ODtnODZOY2soUUJdKGYqKik9bzI6TClbNi8wYUdoZXNtM0NPS2ZaaTF2YFc7Yz1BdVU1bkkncDpUTkFLIXpiNClINTVHWUkheihRVVh3XlR4M2FFIilSR0tqQUh2IlxLXUo/STdaMVpiMklMZyZbKFt0bzxUIzR3Tj5rRG8mZSZHRlRJcXohUmJIVWhmJClRSDxRdjtIIXk5MCgpW3Z1Kzx2PDdQdihcKDMjKXleMWoqKio9KXBSXXBMJT0hby1nKXlwZXFuQTJbISpwal1pLHRCI3kyTz1SXC0hW2JAOkgmZidRYjlZZyFwXHUpNFRsYlVCJj13TlZLb3hORHcjb0BmeSkpUngnWypmZWBJMVIxVXhUcHI2UFhdYi94YyFwLFtLNSI+QGcoND5eaWU7OUdDNmU0JjNUNiYqNDInei4/J2Z0MT9KJUdqNz0ieilbJFFZdHEieiYzQy4pKnpnXT9JbV5zV0x3UHMscHl4bjFGLHhGIWVrYioqXDNWXiY0JFsmeU00J1JQQHM4d0o+bWopMyRScExnJSo9WFN2K10lWyNmXCI+TTxSKXlaIjN1IypwXm1ZMm9BNFk3Xi1gX11EQCEqSFkiW14tLEQoW001KiopM1lfVSJ6blJ6ZWcrai0uIlEnKT54RXEnPi0zIkgiKlJJbldlSjEmKXpFP2Nfblo/RVRUJj5qWDEpKkhcKVw1I0dVYjZcRStUbypRJ3lgNkU2OCgzYUQncCFbL0EpKkgpcD8lWz5vRSkqW1JfLDIoZT9RbVA/bjE5am5sQCM0SGJSRF9TJXBCdl0qNFM7N1cpKikpPS4hUlRjYm8kUTVoQCM+bTMmUitTOkthNno2OCJwZGNNWCxpSFw5WjpGQD1UOCJSNU8ydktqVlp0PTUoSFMjKioqUjQjKXlBPXdZLkxrWnU0ITNdelwieUhFZSM+Zj4lNF53OCRca1hkNThGPzVNQDtOXmNhWTQmbyx6NW1NK1NVZykqPltoTigqUWNHMzYhNEAiKT5mcjxtPiI0WF4lSChRPm1XSSFbTitEW2o8Ink4eGpdRFBjMXApb3kmZitTKGU/S3BFVkFvby88cE9nVTV6VWpzI3ovIipSQjdvdChcc3U/MC51VyZSKzUjZVNKeSs9UkQmUmNadTpQNyElKz8=

evalf[33](HeunD(a, b, c, d, 1/3));

JCJCQGVTSnkrPVJEJlJjWnU6UDchI0s=

myHeunD(-0.9, 9);

JCIrPmkkZXAjISIq

evalf(HeunD(a, b, c, d, -0.9));

JCIram4kZXAjISIq

evalf(HeunD(a, b, c, d, -0.99));

Warning, breaking after 2000 terms, the series is not converging

SSp1bmRlZmluZWRHJSpwcm90ZWN0ZWRH

myHeunD(-0.99);

JCIsIUcmZXZuJSEjNQ==

myHeunD(-0.99, 500);

JCJgamwiKSl5SShHKFItZy9tJSllS3o7VG1cX2glPU4qPUtJYydwaXg/dzc/VSYqKTQlW2VsaEIhcCpcVE0sZS08WyhlZFhuJTM8TjwiM2xObDBnSD5hLzNDKGYhMy0xTXhkKSkqR2V4QiFIJ1FobiRIM04ib0VaPlVCYGYhW015MSEpUlJsdGlCcmVmJUhtOCVcVjVeJik9QSkqKilbZyd6b0l3YydcJno7cmh2JUdeQUUlUmwvKCk0MHZmUHppQV8nNCk0XlRwaCMzbEVxKVIjPjZHSSE9V24vKCpRamxpOkNfT05AKSpIZixUIyp5Mlp1NUNCPCRvaEBjKFJcI2UkPVA/Iyo0L2NCLmxsMDg5a2hyamsiWyEqbyd6X2V2biUhJCsm

random_diffeq := proc(ord, d) uses RandomTools; local myrat; myrat := 'rational'('denominator'=60); { add(Generate('polynom'(myrat, z, 'degree'=d)) * diff(y(z), [z$k]), k=0..ord), seq((D@@k)(y)(0) = Generate(myrat), k=0..ord-1) }; end proc: diffeq := random_diffeq(3, 2);

PCYsKiomLCgjIiM4IiNJIiIiSSJ6RzYiIyIiKSIjOiokRioiIiMjIiIoRihGKS1JInlHRis2I0YqRilGKSomLCgjISIqIiM/RilGKiMiI0hGKEYvIyEiIiIjN0YpLUklZGlmZkclKnByb3RlY3RlZEc2JEYzRipGKUYpKiYsKCMhI1YiI2dGKUYqIyIjXEZIRi8jIiM2RihGKS1GQTYkRkBGKkYpRikqJiwoIyEiKEY/RilGKiMiIzxGKEYvIyEiJCIiJkYpLUZBNiRGTUYqRilGKS8tRjQ2IyIiIUZnbi8tLUkiREc2JEZCSShfc3lzbGliR0YrNiNGNEZmbkYxLy0tLUkjQEBHRlxvNiRGW29GMEZeb0ZmbkZG

evaldiffeq(diffeq, y(z), (1+I)/5, 40);

XiQkIkhqJ1s6SmZaIjNMKj04Vnl3KFtkYlslISNTJCJII0cpeSo0ImVzKCp5LzBgIUdASSkqWyEqPkYl

Displaying the result using the Maple GUI is slow %evaldiffeq(diffeq, y(z), 1/5, 1000000); deq := (1+z^2)*diff(y(z),z,z) + 2*z*diff(y(z),z);

LCYqJiwmIiIiRiUqJEkiekc2IiIiI0YlRiUtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUYrNiQtSSJ5R0YoNiNGJ0YnRidGJUYlKiZGLkYlRidGJUYp

evaldiffeq(deq, y(z), 1/2, 20);

LCYmSSNfQ0c2IjYjIiIhJCI2KysrKysrKysrKyIhIz8mRiQ2IyIiIiQiNUA7aCEzKzR3a2olRio=

transition_matrix(deq, y(z), 1/2, 20);

LUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMvSSQlaWRHRiciKW8mb24j

evaldiffeq(deq, y(z), 5/4*(1+I), 20);

Error, (in rewrite_path) evaluation point outside the disk of convergence of the differential equation (try specifying an analytic continuation path such as [0, 5/4+(5/4)*I])

evaldiffeq(deq, y(z), [0,5/4*(1+I)], 20);

LCYmSSNfQ0c2IjYjIiIhJCI2KysrKysrKysrKyIhIz8qJl4kJCI2JXo7ZCY9Yj5YdzgiRiokIjVfdWlrQSFSY0xeJEYqIiIiJkYkNiNGMUYxRjE=

transition_matrix(deq, y(z), [0, I+1, 2*I, I-1, 0], 20);

LUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMvSSQlaWRHRiciKUdWUWQ=

infolevel[gfun] := 2: evaldiffeq(deq, y(z), [0,11/10], 20);

step_transition_matrix: 1/2 --> 11/10 (1 derivative[s]), prec~=.33333e-21, #terms=96 step_transition_matrix: 0 --> 1/2 (2 derivative[s]), prec~=.30264e-21, #terms=86 absolute_precision_warning: Recall that gfun:-NumGfun works with *absolute* error.

LCYmSSNfQ0c2IjYjIiIhJCI2KysrKysrKysrKyIhIz8mRiQ2IyIiIiQiNVUwPFZ1bUUiKUgkKUYq

TTdSMApJNVJUQUJMRV9TQVZFLzI2NzY4NTY4WCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIjYrKysrKysrKysrIiEjPyQiIiFGKyQiNUA7aCEzKzR3a2olRikkIjUrKysrKysrKyshKUYpRiY=TTdSMApJNVJUQUJMRV9TQVZFLzU3Mzg0MzI4WCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIjYrKysrKysrKysrIiEjPyQiIiFGKyQiNllRS3oqZWBFZlRKRilGJ0Ym